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SECOND-ORDER METHODS IN INVISCID 
SUPERSONIC THEORY 
By J. FELL and D. C. M. LESLIE 


(Sir W. G. Armstrong Whitworth Aircraft Ltd., Coventry: 
Armaments Division) 


Received 21 October 1954] 


SUMMARY 
The second-order perturbation equation for the velocity potential is written down 
for the steady supersonic flow of an inviscid polytropic gas. A particular integral of 
this equation is obtained for the general three-dimensional case, and it is shown that 


this agrees with the known result of Van Dyke for the axisymmetric problem. 


1. Introduction 

In a series of papers (1, 2, 3, 4) Van Dyke has studied the application of 
second-order methods to the solution of problems in the steady supersonic 
motion of an inviscid polytropic gas past one or more slender bodies. From 
this work it is clear that this approach gives results of very satisfactory 
accuracy for many common wing and body shapes, but so far practical 
methods of calculation have only been available for two-dimensional wings 
and for bodies of revolution at zero yaw. The work described in this paper 
makes it possible to develop second-order solutions for a much wider class 
of problems. 


2. The perturbation equations 
We shall briefly deduce the relevant second-order equations, referring 
the reader to Van Dyke (1, 2) for a detailed discussion of their validity. 


Under the assumption of steady irrotational flow, there exists a velocity 
potential Q(x, y,2) which satisfies 


(c*d QO.) 
Bo; Se; ij ; F 1) 
a2 2 (y—1)(2, Q, s 


We use a summation convention and define 


co) 


y is the adiabatic exponent, 
C,. is the velocity of sound in the free stream, 
U is the velocity of the free stream. 

We shall define the quantities M and B 


M U Ic... B 


{Quart. Journ. Mech. and Applied Math., Vol. VIII, Pt. 3 (1955)] 


5092.31 
2.31 Ss 
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D. C. M. LESLIE 
If we approximate to Q by the expression 
Q. = U(x+ d) (2.9) 


and assume that ¢ is so small that we may neglect all terms which are 0 (42), 
(2.1) reduces to the linearized equation 


L\(¢) Pyy T $.2—PB*bex = U, (2.3) 


where suffixes denote partial derivatives. This is not sufficiently accurate 
for many purposes, and so we consider the second-order approximation in 
which we write 
QO=l (x d t ). (2.4) 
In this expression we assume that ¢ still satisfies (2.3) but that we now 
neglect only terms which are 0(9?). Substituting in (2.1), we find 
L(¢) M*{3(y—1)M?(2¢,- $; bi)b;; + 24; $2; +4; $; %;;|, 25 
where i and j are dummy suffixes and denote derivatives. 
A careful examination of orders of magnitude shows that we may drop 
from the right-hand side all triple products involving streamwise deriva 
tives: thus we may write 


L() = 2M*[(N—1)8°4,4¢,,+¢,¢ 





cy | $. br: I a($;, dyy4 $2 $.:) T d, $. 9y2), 


Y2J 


(2.6 


where ¥ =. Yt _.. 
2p 





This equation for ¢ is linear, and, if we can find a particular integral (p.i.), 
the complementary function will satisfy the linearized equation (2.3) and 
can be treated by standard methods. 

Before attempting to find such a particular integral we shall make the 
further simplification of neglecting the triple terms. This is the correct 
procedure for planar systems, and in many other problems the use of an 
approximate particular integral for the triple terms can be justified. 

We finally reduce our problem to that of finding a p.i. of the equation 


L(¢) 2M?| (V al )B*d r Pex 1 py Pry I d. dp re |, (2.7 
where L(¢) = 0. 


Stated in this form the problem is independent of the derivation sketched 
above. We may write the general solution of 


L(¢) = 0 


, f(%, Y.7) 
R—R, 


; ia 
in the form d= —- | dr, 


> 
“4 

. 
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where we are using the hyperbolic vector notation 
a0 A.B 1, B,—B(A, B,+A, B,), 
A /(A.A), 
l f(x,y, un arbitrary function defined over the region %,. 
; [his does not restrict our analysis to those cases in which ¢ is generated 





ysimple source distributions, since, if we allow /(x, y, z) to include singular 


> nctions such as the Dirac delta function and its derivatives, the field due 
=e » system of transverse doublets, ete., can be included in the same 
ia ulation 
[Thus (2.7) can be formally written as 
| He) = [far fff drag eRe 29) 
. “|R—R,|{R—R, 


provided that we can find a p.i. A(R, R,, R,) satisfying 


l 
L(A(R, R,,R,)) ree er ee (2.10) 
R—R, R—R, 
[ W p.1 ol 2 9 S 
||| dv, || | dr, F(R,,R,)A(R, R,,R,). (2.99) 
3, Evaluation of the particular integral 
Since both sides of (2.10) are invariant to a shift of origin, we can write 


A(R, R,, R, u(R—R,, R—R,). 


Furthermore, since (2.10) is Lorentz invariant and the right-hand side 


his equation transforms as a tensor of order zero, we can further reduce 


+ 
S to 


\(R, R,,R,) = »(|R—R,|, |R—R,|, |R,—R,}). 


(3.1) 


his relation has been used to simplify some of the algebra of this section. 


1] 
nee 
In two dimensions, a p.i. can be found without difficulty since the intro- 
tion of characteristic coordinates allows us to write the analogue of 
9) ac 


F(A, [-), 


tecrated directh 
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The analogue of this procedure in three dimensions is to use Whittaker’: | where 


general solution of the wave equation (5). Let 
J(@) 


, = | f(w—By cos 6—fe sin 0, ) dO 


d, = g(x—By cos w—fzsinw,w) dw is finite 








—# if we ] 
where f and g are arbitrary functions, be two solutions of the linearized | '*°F* 
equation (2.3). It may then be readily verified that the p.i. of Inte; 
. . Ie © i 
is — (2) = d1¢2 
sili | “db | des cosec? \(0—w) F(x —By cos 0—Bz sin 6,6) 

x G(a—By cosw—Bzsinw,w), (3.3) | where 
4 
where F(r,s) = f(t, s) dt, Inteor 


and similarly for G(r,s). Now Jones (6) has shown that 
] % dé | 
27 | x—By cos 0—fzsin 0 


Zz. where 





- a —— - 
° if B(y’+-2°)* > a B(y"+-2*) Defini 
] = . 
seso:~C @ < —B ly? +z*)!. (3.4 
Sa2 —B?(y?- 2*)} 
From this result we can evaluate f and g in our special case, and can write; we ha 
l c ¥ >(O—w , , 7 os : 
ene | dw cosec?{| ——]|log(X,—BY, cos 6—fZ, sin @) x 
87128? 2 : 
“7 om But 
< log(X,—BY,cosw—BZ,sinw), (3.5 
where X, = x—2,, etc. This result is only formally correct since the integral 6° the 
has a ssi anual ese’ at 6 = w. However, the expression \ 
»>[9—w : : ee 
x _ } dé [ deo cosec” : Jlog(A 1 — BY, cos 6—BZ, sin @) » 
; = 5 
x log(X,—Y, cos @—BZ, sin! 
3 : a Maki 
cancels the singularity and is formally a complementary function. [hus thr 
; ‘ : ee 


the expression 


A= =| log(X,—fY, cos 0—BZ, sin @).J(9) d0, 
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J(¢) = cosec* ——9 log(X, BY, cos os BZ, sin w)— 

log(X,—fY, cos @—BZ,sin 6)}dw, (3.6) 

finite and is a p.i. of (2.10). The evaluation of this integral is simplified 

fwe put Z, Z, = 0 and use the group properties of the solution to 





store the general values of Z, and Z, in the final answer. 


| Integrating by parts 
] } * cot ( w)sin w do» - Oe 1+ “ 
a As 3Y, COS w “as+t* 
X ye 
az  . t tan 1@. 
' \ BY. i 
oTa parts iva 
l - + 2t 2t 
tan-'t—tan~! \ : ———- | dt, 
27p* . a.J\1+-t® ajz+t 
, X,—BY, 
' X,1+ BY, 
l | 
ng 1st Jedt, 
1 +-t? az+t, 
1 
| _— 
3 | 9'(s) ds. 
: 
| | | | | t* dt n ay ] 
\] 2/1 +s lt+a,s 1+8, 
\ ) ro 9 >y»1 ] > 
oP. ) Xi—BPYT}{X3—B2Y3}# | —— log 2 
he - 
9 ) l , i > 9 ) 
o | X X?— py?) — log|X,+{X5—p*Y3} 
2 26" 7 * 
Making t bvilous gene! tion to Z, Z, ~ 9, we see that the last 
ter! plementary functions while the first becomes 


og (R—R,).(R—R,)+|R—R,|[R—R,|| =A, (3.7) 


m 
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Method 2 
' lf dt 
From the lemma . = 
at?_b (ab) 


we can write (2.10) as 














l ] . dt 
L(¢) — | “Ci © ; {0.5 
R—R,|/R—R, a J (1+#)[(R—R,)?+ 4?] 
(2-4 t 
where R, ~ ——___= i, = — |R,—R, 
1+? 1+? ‘ 
But, using the fact that the solution of 
/ l , a 2 u 
L(x) = REAP where R J(|R/?), 
is a function of R only, we readily derive 
1/A R R24 A? 
tan—“— - 41697 —__— l (3.9 
x al ate 
and therefore 
1 fil » ,) di 
A _ — ti 1Q+1log(1+Q?)—1 ; 3.10 
ap | ig Steet ee o 
— 242 _1 2 _|_ 42]___ 42 2)\1 
where Q = iL(R A " R)*IL Ba ee, =~ (2, 
| t?(R,—R.,)? 
| - l dt 
Define I (8) = Re | ot" 1s +} log(1+-s?Q?) st; r 
, ' 
and then A= | f'(s) ds. 
0 
But 
f'( ) S( | S) 7 Q? dt ] S ( B? ] )? | 
s . i isi iain Se, nes eee enien, Rot MS A 
é ar fi2 | (ij s?Q?)(1 1 ¢2) B?| 1 - (1 s)( B28? | 1)* 
in {| R—R,|+ |R—R,|}*— |R,—R, |? 
where Bb? Sesto | <1. z 
R,—R,|? 
Hence 
1 
' ee 
A f'(s) ds Soe R Ry|+ R R, | : R,—R, A, (3.11 
© : B? ; 2 R, =z. 
0 


neglecting an irrelevant constant. The connexion between (3.11) and 
(3.7) will be explained in the next section. 
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4, The looping technique 





One further difficulty still remains. We have formally written our p.i. 
||| dz, || | dv, F(R,,R,)A(R, R,,R,), (4.1) 
re Ais defined in (2.10). If 6? is negative, then the integrations in 
2.11) are carried out over the whole of ©, and &,. In this case, the limits 
ntegration are independent of R and formal differentiation under the 
gral sign is sufficient to show that (2.11) is a p.i. of (2.9). 
supersonic problems, however, 8 is positive and X, and X, are typically 
surfaces in 0. The obvious procedure in this case is to 
the r¢ yn of integratior to those points on &, for which (R R,)? > U, 
1 simil for X&y, but, since the limits of integration now depend on R, 
lifferentiation gives contributions from the limits, and, in fact, these 
We ca d these difficulties very simply, and in order to demonstrate 
met! ve shall assume that #(R,,R,) is a single-valued, analytic 
nection of ind a, and that &, and &, are the half-planes z, = z, = 0, 
() 0. These restrictions are not essential but will serve for pur- 
ses Of ll tration 
we would expect to take the x, integral from 2, 0 to the 
= v—B| (y—y,)?+-27}}, 
isimilarly for 2. If we could replace these variable limits of integration 
he present difficulty would disappear. This can be done by 
garding nd x, complex variables; a suitable path of integration 
S then a fixed contour C, which loops the branch point and 
passes through the origin at 2, 0, 
Since the path of integration is now independent of R, the operator ZL, 
ied gives 
|| as, | | as, F(R,R2) (4.2) 
LF Lp R—R, R—R, 
where the integrals are defined over these loop paths. 
Since /(R,,R,) is an analytic function of x, and x, and |R—R,! merely 
nges sign on looping ranch point, (4.2) is four times the real integral 
er the re ns of integration defined by 
x, (R—R,)? 0, 


+, (R—R,)? > 0. 


he paths C; and C, are, however, the fundamental ones, and it is only 
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when the singularity at the branch point is of a sufficiently low order tha 
we can use real paths in (4.2). The looping technique can also be used whe; 
we have non-integrable singularities at the branch point, and it is the 
equivalent to the method of finite parts devised by Hadamard (7). 
Thus, if (4.1) is interpreted on a loop path, it is, at least formally, a pi, 
of (4.2). Introduction of these loop paths shows that the four functions 


R,|+ R—R, i R, R,; 








A> pa —_—- -log i | 
: B 2|R,—R, \ 

4.3 
1, ||R—R,|—|R—R,|+/R,—R, ne 
Ag, Aq - log |————— BiB. FA et een 
oe 2R,—R, | 


are all p.i.s of (2.10). We note that d,, is given by 


r 4(A,+A,). 


m 


If we use the special form 
A. L(A, rs T As | A,), (4.4 


then, on looping the branch points, A, merely changes sign in the same way 


as 


|R—R,||R—R,| 


Thus if we use the special form A, in (4.1), that integral too may be taken 


along real paths. A, also has the useful property that it vanishes on the 
two Mach cones 


R—R,| = 0, R—R,| = 0. 


One further point is worth making. Since L is a linear operator, we may 





differentiate both sides of (2.10) with respect to y,, say, and obtain 
{ cA ) B*(y—Yy ) 
“\ay,| R—R,)R—R,|’ | 

We may similarly integrate both sides of this equation with respect to 2, 


over a loop path. Using such techniques it is a simple matter to extend 


our p.i. so that we can deal directly with those cases in which the field is 
generated by multipoles rather than by simple sources. 


5. Relation to known solutions 


It may readily be shown that 


L( M*¢4d, = 2 Ms - Bd - D py { ¢, db ry db, db, : : ; 
which allows us to reduce (2.7) to the simpler equation 
L(¢) = 2M*8?N¢,. ¢,,- (5.1) 


Van Dyke has obtained a particular integral of this equation in the ax! 


symmetric case, and we shall show how this may be obtained by our method. 
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In axisymmetric flow the potential is given by 
| f(&) dé 
} (@—2P—PP} 
Uw2N ae , 
so that < 1 f'(E)f"(m)Ala—€, x—y, | dédn. 
crating by parts, and using the symmetry of the integrand in € and », 
ive 
j WN , CA CA, 
) - 1) SOS (n+) dédn, 
4 ‘ c& c 7) 
he written as 
; “VEE r "'(n) da 
EN | et | (ee 
‘ co p*r*;* J We—-gP ee 
M2A 
| is \ Dyke’s rest 
We h verified our p.i. gives the known result for two- 
mensior und for yawed cones. This will be proved in subsequent 
In principle, the triple terms in equation (2.6) can be treated by these 
ethods and an extension to higher dimensions is also possible. Preliminary 
nination, however, suggests that the functions involved are very 
ymplex 
[he forn e obtained here have been applied to the study of second 
| i] 1 1 1 . . . . 
er flow d yawed bodies of revolution and to the calculation of 
ickness corrections to the aerodynamic derivatives of wings of finite 
spect 1 This work will be published later. 
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ROTATIONAL FLOW THROUGH CASCADES 


PART I. THE COMPONENTS OF VORTICITY 


By W. KR. HAWTHORNE 


(Department of Engineering, University of Cambridge) 
[Received 25 June 1954] 


SUMMARY 
The flow of a fluid with a non-uniform velocity passing through a cascad 
turbine or compressor blades is examined theoretically to determine the natur 
of the components of vorticity in the direction of flow. Three such components at 
found, two lying in the stagnation streamline leaving the blade, and one distribut. 
in the stream. The distributed secondary circulation is that due to the curving 
the flow in a bend, or passage between the blades, as described by Squire an 


Winter. The trailing shed circulation is due to the change in circulation along tl 





blade as found in the flow of a uniform fluid about a wing of finite span. The thir 
component—the trailing filament circulation—is due to the stretching of the vortex 
filaments carried with the flow between the upper and lower stagnation streamlines 
in the wake of each blade. The application of these results to cascades and isolate 
aerofoils is discussed briefly. 


1. Introduction 

WHEN a fluid with a non-uniform total pressure flows around a bend 
secondary flows may be generated which can materially alter the character 
of the flow. Such flows have been studied by several authors (1-5) and 
a satisfactory theory obtained for the flow of a non-uniform inviscid in- 
compressible fluid when the secondary flows are not large. 

These flows are of particular interest in cascades of turbine and con- 
pressor blades where non-uniform approaching flows of considerable extent 
may be found on the walls bounding the cascade. Recently Preston (6), 
using a simplified approach, has discussed the generation of the component 
of vorticity in the direction of flow, which gives rise to the secondary 


flows. Preston starts with the concept that the vortex filaments carried 


with the rotational flow may be followed in the Lagrangian manner 4s | 


the flow passes through a cascade. Kronauer (3) and Hawthorne (5) have 
also made some use of this approach. When studying the flow past 





cascade, Preston assumes that in addition to the components of vorticity | 


present in the vortex filaments, a change in circulation about the blade 
or aerofoil will cause circulation to be shed off into the stream. Th 
existence of this shed vorticity has been demonstrated in the theory 

finite wings in a uniform flow. The concepts are, however, less familial 


{Quart. Journ. Mech. and Applied Math., Vol. VIII, Pt. 3 (1955)] 
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non-uniform flow and there is a need for an analytical proof of the 


existence of the different components of vorticity. Such a proof is 
ttempted this paper. 





Two analytical approaches are made. One is based on the equations 
r the flow of an inviscid, incompressible fluid in a force field. In the 
ther the vortex filaments are followed as the fluid passes through the 
| cascade. By demonstrating the equivalence of the two results various 
he flow are explained. 
) 
2. Secondary vorticity in a flow with a force field 
\n expression for the change of the component of vorticity in the direc- 
n of flo lerived in (1) for inviscid, incompressible flow in the absence 


f body f \ similar method may be used for the case when body 


rticitvy Q curl V, where V is the velocity vector (scalar q), 


vy be split into two components, namely a component 
@.V)., f\ .. 
dy (=) 
\ V) q 
ng a sti miine. where ; the scalar of the component of vorticity in 
how, and a component 
(V Q) V 
(V.V) 
¢ umling Since 
div Q (), 
; . VA(VAL 
t follows that div (EV /q)—div - ae 0; (1) 
(V.V) 
0 nd since div (EV /q €q) div V+-V.grad(é/q) 
: div V 0, 
~. VA(VAR 
Wve V .grad(é/q div — i : (2) 
? re 
mn l equatiol 
a : | ‘ " 
V .grad)\ grad p+-F (3) 
p 
here F is the vector representing the body forces, and using the identity 
orad(V.\ 2(V .crad)V+2V A curl V 


- n (V.grad)V = grad 3q7—V A Q, (4) 
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we have by subtracting (3) from (4) 


VAQ = grad p,—F, 


> |= 


where p, is the total pressure (= p-+ pq"). 
Substituting the expression (5) for V A & in equation (2), we find that 


V .gradié/e) = div! +. (Vv / -eracl ps) _diviiwar\. ¢ 
\¢? p \g° 
Now the circulation round any stream tube of cross-sectional area dA js 
EdA and, since gdA, the volume flow along the stream tube, is constant, 
€/q is the secondary circulation per unit volume flowing. The term on 
the left-hand side of (6) is then proportional to the gradient of the 
secondary circulation in the direction of flow. 


If no component of force exists in the direction of flow (inviscid flow), 


V.F = 0, (7) 


and, therefore, V.grad(”*) = 0. (8) 
p 

Condition (8) expresses the fact that the total pressure remains unchanged 
along a streamline. The second term in equation (6) then represents the 
divergence of a quantity f/qg, where f is the scalar of F. This quantity 
is analogous to the circulation about an aerofoil which the force field may 
be supposed to represent. Applying equation (6) to the case of a wing 
in a flow with p, everywhere constant, the right-hand side becomes equal 
to the rate of change of circulation along the span of the wing (since the 
first term disappears) and the well-known result that the strength of the 
trailing vortex sheet is equal to this rate of change of circulation 
follows at once. 

The first term on the right-hand side of equation (6) may be expanded 


as follows: 


\ 


div! a(V grad (9})) = 3. grad (”2) Vlerad (“,) grad (2)), (Sa 
\q° p//} @ p \ 7° p}} 


Substituting for grad q? from (4), (8 a) becomes 


] . 
- Q. grad (2 . 
q" p 


I ae Ps)) . oe y l Loar ; _— Po | u Q (Rh 
aa i grad | tab .grad)V - 2V grad (7). 
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Expanding the last term and noting (8), (8 b) becomes 





42 grad (72 I foy grad (P*) .(V.grad)V) 


qi p 
| 7 | 2(V-V)(graa (22) .2)} (9) 
q' | p 
The last term may be rewritten, using (5), in the form 
| srad (2°). — 1 (VA, 2.24F.2) = +(F.Q). (10) 
p ] 
Equation (¢ iv therefore be written 
V.orad | LI 2 srad (£2)! .(V.grad)V 'F.Q div (* nat 
q* | p/) q° % og 
(11) 
In the absence of body forces this equation reduces to that derived 
ously by the author (1). It describes the generation of a component 
ity in the direction of flow in a curved duct, first noted by Squire 
| Winter (2) and later discussed by other authors (1, 3, 4, 5). Ina 
scade the first term on the right-hand side of equation (11) describes the 
| component of second uy circulation due to the bending of the flow in the 
ce between the blades. From this must be subtracted a term (the 
he right-hand side of (11)) which represents the rate of change 
nd the blades in their spanwise direction. This may be 
cht of tex sheet trailing downstream from the blades similar 
that obt | with wings of finite span in uniform flow (p,) = constant). 
In non-unif flow, however, another term (the second term on the 
right-hand side of (11)) must also be subtracted. The physical significance 
this te | be shown later 
3. Application to cascades of aerofoils 
Conside nitially quasi-two-dimensional flow (Fig. 1), in which the 
ssure nitially uniform and the initial velocity U, varies in the 
lirection y. Let the components of the velocity vector V be C,, C,, 
LC, in j-, 2-directions respectively. Let this flow approach and 
deflected |} . force field (e.g. a two-dimensional cascade) of finite 
ith 1 lirectior | of infinite extent in the y-direction. In the 
lirection similarity will be preserved so that all properties of the flow 
| lep ly 
Suppose, first, that we envisage a flow in which the original planes of 


nstan emain planes as they pass through the force field, i.e. pp» is 
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R. HAWTHORNE 
a function of z only. Also, let the angle « made by V with the z-axis }p 
independent of z. By continuity there will be no acceleration in the +. 


I 


direction. Hence C’. will be independent of x and y and «a = tan-(C_/¢ 
“ull, 





Cz 
a 
A] 
Ss } 
a nth —-— x Cx | 























will be independent of y and z. The secondary circulation may be obtained 


directly from equation (2). Here 
V =i€,4+jC, = ic,+jC, tana, (12 
: . a, . .dC, dtan 
Q = curl \ —itana—“74 j—*+kC, “i (13 
dz dz ; dx 





where i, j, and k are unit vectors in the 2-, y-, z-directions respectively 
Also 
' a dtan . 9 dtan = GG. 
VAQ = iC? tana ; {Ce - a“ kC,,sec?a —, (14 
dx oC eee dz 
. a — ene} . d P 
VA(VA&) = iC? tan a sec?a —* + jC? sec?a —2 —kC3 sec*a- =, (15 
, ill dz dx 
and 
VA(VAQ VA(VAQ : dC. : dC. , da : 
bas Eh, a = itana—2+jC,—_* —kC;,sec*a—. (6 
q* C2 sec?a dz ~ dz dx 
Hence 
~ (VA(VAQ) dc, > da dC, ., da Je 
ihe —_— | = — S€C“a — — —— sec*a Q. (17) 
q- dz dx dz dx 


Hence the component of vorticity in the direction of flow which is initially 
zero remains zero in this type of flow. 
Such a force field does not, however, correctly represent a two-dimen- 


sional cascade since, in such a cascade, there is no component of foret 


exerted by the blades in the z-direction, i.e. along the span. There wl 
therefore be an acceleration in this direction, since in any surface © 
constant py 


P 5. 2.9 PX | 1 (18 
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Furthermore, g in any such surface tends to be proportional to U,, and 
there will therefore usually be a component of the pressure gradient normal 
to the surface. In the absence of any other force this component of the 
ressure gradient will cause an acceleration and the surface of constant p, 


no longer remain plane but will be displaced in the z-direction. Such 





splacements have been observed in the flow through cascades imme- 
tely upstream and downstream of the blades as well as within the blade 
passages 
In order to develop a first-order theory it will be assumed as in references 
2) and (5) that the basic flow at each value of z has a constant axial 
itv C. and an angle independent of y and z. The velocity com- 
ents of this basic flow are then given by equation (12). Superimposed 
such a basic flow will be small perturbation velocities which are so 
ll that the value of p, at each value of z remains constant. This sup- 
sition affects the value of the vorticity component, 2,,, the com- 
nent of Q in the xy-plane and normal to the velocity vector. In the 
w described above, with no component of acceleration in the z-direction, 


) 


3) and (5) is given by 


Irom equations 
gQ 4. Ms. af. 


(19 
, p dz p dz 


ich agrees with the value deduced from equation (13). 

In the absence of a force in the z-direction equation (5) yields the result 
qQ dpo 
a «zy 


p dz 
‘ 


U, du, 
q ‘dz . 


(20) 


is is different from the value deduced in equation (13) from the com- 
nents of the basic flow only. This difference may be explained by show- 
g, as in reference (5), that the component of vorticity due to the small 
turbations is of the same order as the vorticity in the approaching 
w and may not, therefore, be neglected. 
ipproximation the force f acts in the xy-plane normal to 
the basic two-dimensional flow and is approximately 


he vorticity component in the z-direction, by the equation 
f = qQ,. 21) 


et Che force exerted by a cascade of blades on a non-uniform flow is not 


eadily predictable. The flow presumably leaves the trailing edges of the 


LdeS SO 


that.] 


ukowski’s hypothesisis satisfied, or approximately satisfied. 
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Hence this condition becomes one of the boundary conditions of the 
problem and enables the blade force to be determined. 

The vorticity component 2, comprises a component due to the basi 
two-dimensional flow and a component due to the small perturbation 
which may be of the same order. A knowledge of the perturbation velocitie: 
which satisfy all the boundary conditions is therefore needed to calculate () 
and thence f. As the perturbations may be regarded as due to the additiy, 
effect of the secondary components of vorticity generated in all the surfaces 
of constant total pressure, the effect of one surface on another must } 
considered. This problem is considered further in Part IT of this paper, 


For the purposes of this part of the paper it will be assumed in bot! 


force field exerted is that corresponding to the compenent of Q, due t 
the basic two-dimensional flow only. This implies that at a given val 
of z the circulation about each blade of the cascade will be that whic 
would be obtained in a two-dimensional flow with an approach velocity 


valid in that it is possible to conceive of a cascade of twisted blades 
producing such a field of force. 

Alternatively, a limiting cascade may be envisaged in which closel 
spaced blades fix the outlet angle of the perturbed flow at a value ink 
pendent of z. The perturbation velocities are, thereby, controlled in su 
a way that the component of vorticity, €, at outlet from the cascade for 
field can have no value other than zero. In other words, one surface | 
constant total pressure cannot slide over the other. The force field require 
to produce this effect will now be different from the previous one. This 
example is, however, not as useful for analysis here as the one chose 
above. 

The acceleration of the basic flow in the xy-plane is C®(d tan «/dz) in 
the y-direction, whence the force is C2 sec a(dtana/dx). Hence, inserting 


appropriate values in equation (11), we find 


‘ d(E/q) 24 U, dU; C2 cos _dtan v l C2 sec _dtan a U, du, 
ds q dz dx a dx dz 
d (C2 seca dtana a 
- dz q dx : 7 


where ds = dxseca is an element of length in the direction of flow. 
Substituting g = C,,seca and noting that « is assumed to be independent 
of z and C,, to be independent of 2, 

rae | da |. yy dU, dx , AC, secada (22 


l dz 1 dz Gg © dz C re 


I r 


d(é/q) 2 7 - 
2 sec?y 





| 
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; ; ; _ 
analytical approaches that in each surface of constant total pressure th: 


equal to the value of U, corresponding to that z. This assumption is | 
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f +] [he first term on the right-hand side of this equation is identical with the 

expression obtained in (5) for the secondary circulation in a similar flow 

bas ) it am ibsence of body forces. It is the expression which is required to 

at determine the secondary circulation obtained by the bending of the flow 

cit etween 1 ides of a cascade. The integration of this equation will 
ite() | be performed after the discussion of the vortex filament theory. 


4. Vortex filament theory 


st },? Fig. 2 shows a surface of constant py passing through a section of a two- 
™ ide. It will be assumed that the basic flow in this surface 
4e 





aah 
» od 
ea | U, 
—— | F 
Pi Aa 
Fic. 2 
| | 
limensional potential flow and that flows due to the varia- 
i mal to the paper may be regarded as perturbations of this 
w. The analysis of (6) and the earlier sections of this paper may then be 
| appled to the flow. In the flow there will be vortex filaments moving 
th the fluid and satisfying the theorems of Kelvin and Helmholtz. One 
S unent, ab, lying in the plane of constant po, is shown in 
Fig. 2. The vortex filament will always be composed of the same particles 
| and its motion may be followed until it is downstream of the cascade. 
Vwing to the differing velocities on the convex and concave surfaces of the 
les, 1t will take up a position def, Fig. 2. The filament may stretch 
irom @ to e along the stagnation streamline or, as is suggested in 
6) for 1 led leading edges, the filament may stretch upstream from 
round the nose of the blade and back along the convex surface as is 
i shown tor the lower blade in Fig. 2. The portion ef of the vortex filament 
7 
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contributes a component of vorticity in the direction of flow. Th 
secondary circulation due to this vorticity component in a plane stri 


normal to the flow stretching between streamlines %, and ys, and of unit 


height normal to the paper is, from equation (2.4) of (5),+ 


60, = | (l 1 Vibe dd — (l 1 Win i, (23 


where ¢ is the potential function for the basic two-dimensional motio 


and dd (q/U,) de. 94 


The integral is in each case taken along a streamline. Denoting by g 
the velocity on the stagnation streamline y%, taken along the conv 
surface of the blade and by q’ the velocity along stagnation streamlin 
ys, taken along the concave surface of the blade (Fig. 2), the equatior 
becomes . ” 


SI, aU} [ (3) ag | (=) ag}. 25 
dz\ J \q’ J \q 


¥ x 


Far upstream and downstream the potential lines such as ab and gh wi 
be normal to the flow. Taking the integrals over such limits and subst 
tuting for dd, we obtain 


h { 
P - dU; | [ U, ds r eh ds| 
 da\|j oq’ q | 


a b 


The equation of (5) needs slight modification. It was originally applied to a strut 


isolated aerofoil for which U U,. The modification for U, U, is as follows: 


From (5), equations (1.18), (1.20), and (2.3) vield the downstream vorticity, 
1 


> U2 dU, l c(q/U,) 
2 = ——— ds, 
l 1 {Zz i (q if )” con 
vhere the integral is taken along the streamline and n is a coordinate normal to the stt 
line lying in the plane of constant pp. Now 
uso 
TV, | é, dn, and dd (q U,) ds, dis (q U,) In. 
1 
Hence downstream at infinity 
us 
U, 
‘ ‘ Lil 
a, 7 ; 


On substitution for dn and ds, the expression for €, may be used to obtain the expres 


for ols: equation (23), by integrating once with respect to yw. 
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Because of the periodic repetition of the flow in each blade passage, the 


Raut ¢ 
ist 1 


ntegral has the same value as a similar integral taken from c to e 


long the stagnation streamline passing along the concave surface of the 


7 | 


. 
a a h 


u ds ae al. 


The first two integrals cancel out except at the blade and the last two 
v be evaluated from the geometry of Fig. 2, thus giving 
} a L}] | ds 1 osin Xy a(U, U,)sin mi (26) 
dz\ "J q 
re o is the pitch of the blades and > _ is taken from the leading edge 
J 
ong the convex surface to the trailing edge and back along the concave 
ce. It is therefore a negative quantity. 
Noting that, by continuity, the downstream velocity U, is given by 
LU), cos a, U’, cos Xo, (27) 
aU, | ; “ds osin 2d sin 20, o 
n 75 | é, p = Sanne ——|. (23) 
. a q kd COS Xy 
positive direction of 5, is that of a right-handed screw advancing 
nstreal 
Each vortex filament ha portion de (Fig. 2), which produces a secondary 
tion opposite in sign to that produced by fe. Only one filament is 
wn in Fig. 2 but, in fact succession of them could be drawn between 
nd e (one such is partly shown at d’ and e’), all of which will contribute 
the circulation about the stagnation streamline at e. The circulation 
is the m of the cir tions in all the filaments such as d’e’, where 
es between d and c. The filament starting as ab reaches de in a given 
Since the time ¢ taken by a particle starting at a and reaching d via 
concave surface of the blade is the same as the time taken by a particle 
ing at nd reaching e by the convex surface, we find that 
d ¢ 
, [ * ds [ ds | ds 29) 
q Ja Jaq 
a d 
supe pts ‘and “ refer to the concave and convex surfaces respec- 
ly. Co er a filament such as d’e’ which starts from a’ so that d’ 





e same tims Then the total circulation at e will be the 
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circulation associated with all the vortex filaments between a and a’ and 





since - 
r ds rds 
peal ‘ 
J @ J 
a d 
a’ 
. “ds rds 
and, from (29), | = sd > 
J q J 4 
a 
: _ fds 
or aa U; ) 5 (3] 
J 4 
Hence the circulation about e is given by 
_ , aU, _ au, fds 
or, — aaa - U, l d * . 39 | 
5 dz dz J q 


The equivalence of this term and the F.Q term in equation (11) may b 
approximately established as follows. From Fig. 3 the force exerted by 














an element of blade (assumed to be thin) of length ds on the fluid is given 


by (p’—p")ds 4o(q"*—q'*)ds, since py is constant in the surface. The 

force f per unit mass of fluid is then given by 
f q’ t q’ q’ q’ 

7 ; 


a n 


99 
(00 


Noting that 3(q’+-q’') is approximately the average flow velocity q and 

substituting this expression for F and equation (20) for Q in equation (11), 

the contribution to the secondary circulation from the F.Q term is give! 

by . : woe 

: d(&/q) Vu dU, 1 a(q"—q’) 

J hee oe eee, 
ds dz q n 


By continuity nq is constant and equal to U,acosa,. Hence 


U,, o COS a5 — 
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q~ 


Z py [ (q"—q') ds 

U, t dz 

or ol, U; d L | q f ds. {vv 
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The equivalence of this with the expression for 5I,, equation (32), is 
htained by noting that 
*ds aie | ea’ —d gn 3 ; 
| | ds u : 4 ds = f “4 ds. (36) 
{ me / J 44 ‘ q° 


It follows that the physical significance of the second term on the right- 


nd side of equation (11) is that it represents the component of vorticity 














} 
' 
' 
} 
} 
. 
| 
| 
| 
to the stretching of the vortex filaments along the streamline from 
the trailing edge, as shown by de in Fig. 2. 


lo complete the analysis it is evident from equation (11) that the com- 


ponent due to the change of circulation about the blade must be added. 


his is give 


dl’ 
tan a,)) a COS a 7, (tan Yo — tan a). 
dz 


(37) 


1] sin 9 in? 
_ , , {sin 2~,—sin 2 
| . ecetane Fe x, (tan a,—tan .,)\. 
| ? cos Xy a | 


(38) 
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This is identical to the expression obtained as follows by integrating 
equation (22a) from far upstream to a section far downstream. Noting 
that C, is independent of x, we obtain 


9° 9 
5 U, dU, Ll aC, f 
&. U, = = - d (2 cos*a— 1) dx — — : | sec?a da, 
Sie cy Us dz Cc. a 
a 
1 1 
: l dU : . cos a, dU 
—— ee 1/1(sin 2a4,—sin 2a,)|— — : ' (tan «,—tana,). 
C9 2 2 1 2 1 
COS a, COS a, dz cosa, dz = i 
(39 | 
Hence aie 
él €,0 COS Xs 
dU, (sin 2a, sin 2a, | 
o —— | —__—————-_— €08 a, (tan a, —tan a,) 
dz \ 2 COS a, ‘ 
dU,. ‘COSA COS a, 
o = SIN ao — — ——}. (40) 
dz “\cOsa, COS a, 


: oe a) ee : ; 
Values of 4] |, 7 ! from the above expression are shown graphically i 
/ az ; 
Fig. 4. 
5. Isolated aerofoils 
For an isolated aerofoil the vortex filament approach leads to the sam 

values for the trailing shed circulation and the trailing component of the 
vortex filaments as those given in equations (32) and (37). Equation (37 
has, however, to be modified. The circulation K about the blade is giver 
by . : 
‘ kK oU, cos a, (tan ¥,—tan «,), (41 
so that, for the isolated aerofoil, we obtain 


. K dU ' 
si, = —— —-, (42 

U dz 
where U is the velocity approaching the isolated aerofoil. An expressio! 
for the net distributed circulation has been obtained in (5), equatio! 


(2.8), and in terms of the sign convention used here it is 





—_ dl’ _f@ee £ 
Sr; Ses as 
dz J oq l 
The net value of secondary circulation is then 
: sie es , KdadU 
rad | 61,+60,+802, = —— 2 —— ——., (4H 
bi U dz 


i.e. twice the trailing shed circulation. 
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The assumption that the circulation in each surface of constant total 
pressure is that corresponding to two-dimensional flow with the corre- 


sponding value of U must, however, be remembered. It is not valid in 


e theory of wings of finite span and is unlikely to be valid here. 


. Conclusions 


The analysis shows that the components of vorticity in the direction of 
w produced when a non-uniform stream flows through a cascade or 


passes an aerofoil, may be divided into three parts. The first is that due 


the bending of the flow in such a way that the streamlines are no longer 
seodesics on the surfaces of constant total pressure (see (1)). This is the 
secondary culation described by Squire and Winter (2). The stream- 


wise components of vorticity are distributed throughout the flow. This 
ributed secondary circulation is the result of a distortion of the vortex 


uments carried with the stream and arises in flow around bends, through 


the passages between blades in cascades, and in the flow about obstacles. 


The second component—the trailing shed circulation—is that due to the 
hange of circulation about bodies immersed in the flow. This is a familiar 


meept in the theory of wings of finite span immersed in a uniform stream. 


he analysis shows that the concept may be extended to non-uniform 
ws, and confirms the assumption made in Preston’s simplified treat- 
nt (6 

The third component—the trailing filament circulation—arises because, 
passing round the aerofoil or blade, the vortex filament on the convex 


surface moves ahead of that on the coneave surface so that the filament 


to stretch in the wake. The stretched filaments lying in the trailing 
stagnation streamline contribute a vorticity of the same sign as the trailing 
shed circulation and, together with the latter, form a trailing vortex sheet 
eparating the flow over the upper surface from that over the lower surface. 
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PART II. THE CIRCULATION ABOUT THE CASCADE 
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SUMMARY 


The variation in circulation about aerofoils in cascade when the approaching 


stream has a non-uniform velocity is determined by assuming that the secondan 
flow is confined to the rectangular passages between, and extending downstrear 
from, the blades. The strength of the trailing vortex sheet separating each bla 
passage is estimated from the secondary flow. From this is deducted the trailing 
filament circulation, defined in Part I, to give the trailing shed circulation fro 
which the variation in blade circulation can be determined. Measurements 
pressure distribution on a cascade of impulse turbine blades confirm experimenta 
the variation of lift and circulation predicted by the theory. It is shown that t 
secondary flow tends to cause a stalling of the blades near the wall. A sim 


behaviour for isolated aerofoils in non-uniform flow is predicted. 


1. Introduction 

To demonstrate the existence of three constituents to the component of 
vorticity in the direction of flow produced when a non-uniform flow passe: 
through a cascade, a simplifying assumption was made in Part I (1). I 
was assumed that the forces exerted on the fluid at any point on the blad 
span would be the same as those obtained in a two-dimensional flow wit 
the same velocity of approach. While it is possible to conceive of a cascad 
of twisted blades producing such a force, the assumption is clearly not 
generally valid for an untwisted two-dimensional cascade as normall 
used in cascade wind tunnels. Even in non-uniform flow it may be pre 
sumed that the circulation about a well-designed aerofoil in a give! 
cascade will tend to adjust itself so that Joukowski’s hypothesis, that the 
velocity at the trailing edge is finite, will be satisfied. 

Regarding the flow as a basic quasi-two-dimensional flow upon which 
are superimposed perturbations induced by the secondary vorticities, this 
boundary condition at the trailing edges of the blades must be satisfied 
not only by the basic flow but also by the secondary induced velocities. 

The basic flow fulfils the required condition since it is a series of layers 
of two-dimensional flow in each of which the outlet angle and blad 
circulation are determined by the condition that there are no infinitt 
velocities at the trailing edge. The boundary condition for the secondat) 
flow may be obtained by making the simplifying assumption that the shee! 
of fluid leaving the trailing edge of the blade remains, at least for 


‘ 
[Quart. Journ. Mech. and Applied Math., Vol. VIII, Pt. 3 (1955)] 
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ppreciable distance, the boundary separating the flow from adjacent 
wntiae s. Some distance downstream a diagram of the secondary flows 
) ina plane normal to the main flow would appear as in Fig. 1. The 
secondal low between the sheets of fluid trailing from the blades will 
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Bi Secondary flow downstream of cascade. 
similar to that in a passage after a bend and there will be no velocities 
? at the sheets normal to or across them. With the above assumption the 
ondary flow may be calculated by the methods used in (2) for calculat- 
« the velocities induced by the distributed secondary vorticity generated 
the flow of a non-uniform stream around a bend (2, 3, 4). The strength 
f the trailing vortex sheets situated at y= = b, +3), t 5b, ete. (Fig. 1), 
\ in be determined by calculating Aw, the change in spanwise secondary 
elocity cross them 
I : 
ws 
‘ ' 
ALY 
\ y 
liiy £ 
oa 
A = WALL 
U J, ~~ " . 
¢ a / Fs § 
<i _ 
“w 
SECTION AA 
‘1 Fic. 2. Notation. 
8, t As shown in Part I, the sheets of fluid leaving the trailing edges are vortex 
tisi sheets consisting of the trailing shed and trailing filament vorticity. If 
ies the trailing filament vorticity is estimated and deducted from the total 
strength of the vortex sheet, the trailing shed circulation, and thence the 
iriation of circulation along the span of the blades, may be obtained. 
n Chis method has been applied to some experimental results from a turbine 
nd cascade in which such estimates of the variation in circulation can be 
> § mpared with values obtained from measured pressure distributions. 
Che notation used is shown in Fig. 2. Following references (2) and (4), 
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the secondary velocities v and w may be described by a stream functio 


bb (w 6ys/Cy and v = —oys/0z), whence €,,, the vorticity in the directio, 
of flow downstream of the passage between the blades, is given by 
Ou Cus 


; ( 


+ 
cy~ C2" 


Vas 


The value of €,, for a rectangular bend such as an impulse turbine blac 
passage was shown by Squire and Winter (2) to be 

av, 

dz’ i 


where « is the deflexion and U, the velocity far upstream. Equations (| 


:- 9° 
6 > 


ut 


and (2) have been solved for the experimental variation of U, with z to 
give the stream-function % and thence the value of w on either side of the 


trailing vortex sheet (Fig. 1). The solution was obtained by the relaxation 


method. The value of deflexion, «, used was the experimental valu 
obtained from two-dimensional cascade tests in which the approacl 
velocity was uniform. 


2. Theory applied to impulse turbine cascade 
A cascade of (approximately) impulse turbine blades shown in Fig. 3 


of 6-in. chord and 18-in. span, was tested with the non-uniform approac! 








SCALE: reer 
0 6" 12 
Fic. 3. Turbine cascade. Chord 6 in., pitch 6 in., span 18 in. 


velocity shown in Fig. 4. Measurements of secondary velocities were made 
in a plane a short distance downstream of the trailing edge, and pressure 
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distributions on the surface of the blades were measured at various span- 
wise positions. The direction and magnitude of the resultants of the 
pressure forces are shown in Fig. 4. As the flow was symmetrical about 


the midspan position, the curves in Fig. 4 are drawn for only half the span. 




















[he subscript ¢ is used in the figure to signify conditions at the mid- 
span position 
| 
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[he secondary velocities were calculated from equations (1) and (2) where 
| 
he defiexion, was assumed constant at the two-dimensional flow 
ue. Theoretical and experimental values of w, the spanwise velocity, 
compare l at two spanwise positions in Fig. 5. 
The values of Aw, the stre ngth of the trailing vortex sheet, were calcu- 


ited from the theoretical values of w. They are shown in Fig. 6, curve A, 


integrated values, the value at any z being the total strength of the 

wing vortex sheet from the midspan position to the spanwise position 2. 

With the results obtained from the pressure plots it is possible to show 
that the + 7 


ig vortex sheet consists of components due to the trailing 


“ment circulation and the trailing shed circulation. The trailing 
ment cireculat 


ion has been estimated from the theory of Part I and is 
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wn as curve B, Fig. 6. The trailing shed circulation is theoretically 
en bv the difference between curves B and A shown as curve C, Fig. 6. 
is is compared with experimental values of the blade circulation ob- 
ined from pressure plots and shown as experimental points in Fig. 6. 
e method of calculation is as follows. 

The t ¢ filament circulation from equation (32), Part I, is 
’ aU, fds 
Ae p (3) 
=4 dz q 
the distributed secondary circulation from equation (28), Part I, is 
r al, | ] ‘= o Sin 2a,—sin 2a, ) 
, ee q z= COS Xy | 
COS ap. (4) 
Squire a Winter’s (2) estimate for €, given in equation (2) has been 
de and hence from equations (2), (3), and (4) we have 
, il l sin 2a~,—sin 2a - 
: . Ze COS Oy}. (») 
(2 COS a, " 
ue of 60, for values of « and a, appropriate to two-dimensional 
w have been calculated and are shown in Fig. 6, curve B, as integrated 


ies from the midspan position to spanwise position z 

rhe difference between the curves for Aw and 89I, is shown as curve C 

Fig. 6. It should give the total circulation shed off the blade between 

midspan position and z and hence the difference in circulation about 
| the blade at these two positions. For a blade in cascade the circulation 
1 to the lift L by 


y pU,, I, (6) 
Where p is the fluid density and UL’, is the vector mean of the velocities 
fentry and exit, U, and Changes in [, from the centre span value 
have bee termined from equation (6) and the lift and velocity values 
shown in Fig. 4. These values are shown as experimental points in 
tig.6. They are to be compared with the curve C, showing the difference 
etween Aw and oI. Considering the approximations made and the large 
v» and « due to the secondary effects, the agreement is 
\ second set of curves have been drawn in which ¢, has been based, 
n (2), but on the more exact equation (4). The value of 

/ has been determined from measurements of the pressure distribution 


ude force measured from the pressure plots is approxi- 


Wi lade is stalled and the drag is large, this approximation 
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about the blade in two-dimensional flow. The results are shown in Fig, 7 


and indicate a fair agreement between theory and experiment. 
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Fic. 7. Values of trailing circulation (€ from equation (4), $ ds q from 
pressure distribution). A. Trailing vortex sheet from Aw ——. B. Trailing 
filament vorticity o—o—o. C. Shed circulation by difference »—>—. 
Blade circulation from measured pressure distribution x x x. All curves 


show integrated value from midspan to position z. 


The existence of a trailing filament component of the trailing vortey 
sheet is clearly indicated by the results shown in Figs. 6 and 7. In this 
instance it is of the same order of magnitude as the shed circulation and 


is, therefore, about half the strength of the trailing vortex sheet. Experi 


mental confirmation of curve A in Fig. 6 is obtained indirectly by the 


comparison of measured and calculated values of the spanwise velociti 


w shown in Fig. 5. 


3. Modifications to circulation and lift due to secondary flow 
In Part I it was arbitrarily assumed that the circulation about the blad 
at any spanwise position would be that obtained in a two-dimensiona 
flow with approach velocity U,. This gave a value of trailing shed circula 
tion . 
ols —o COS a, al ! (tan o,—tan N4)s 
dz * 
where «, and a, were assumed invariant with z. The circulation based 0 
this equation is compared with the result obtained from the measur 
pressure distributions in Fig. 8. The actual change in circulation is mu 
smaller than that given in equation (7). 
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If the assumptions of Part I were made for the lift, the value of thy 




















a Bs , 4, To 
lift coefficient C,, defined from - 
The 
wo: oe } 
L Spl mL es (8 ABCI 
m : the su 
would be uniform along the span. Values calculated from the pressure 
plots are shown in Fig. 9. The increase of C, towards the wall can ly | 
regarded as due to an upwash. It is accompanied by the increased angle 
of attack a,, (i.e. angle of resultant of pressure forces) reproduced jr 
ft . - nr » we F - wna, 8 
Fig. 9 from Fig. 4. The results of Fig. 9 are replotted in Fig. 10 as 
curve of C, against a,,. Both Figs. 9 and 10 show that very near the wa 
C, continues to increase with decreasing distance from the wall. Th 
L 
angle of the resultant pressure forces, however, after increasing, suddenh 
decreases. 
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The results of tests of the cascade in uniform two-dimensional flow show | any p 
that this profile begins to show increased losses at a C, of about 3» In] 
On stalling it is to be expected that, owing to the increase of drag, the angle that | 
of the resultant pressure forces will become less than the vector mean 0! Howe 
the inlet and outlet air angles normally used to define angle of attack i! the s| 


cascade theory. At large values of angle of attack the reversal of 


sins ; 7 5. Tr 
and OC; in Figs. 9 and 10 respectively can be partially accounted for ! 
this effect. The agreement between the two curves in the unstalled regio! In 
is fair wna 
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©} 4, Total secondary circulation 


The total secondary circulation 6 embraced by a circuit such as 
{BOD, Fig. 1, where AB CD 2b and AD BC = 82, is given by 
e sum of the distributed and trailing vorticities, i.e. 
| €,, 62 dy+Aw.9z, (9) 
since the velocities along BC and AD are the same, by 
B D 
| ? dy v dy 
| r ¢ 
| 
er 
| dz dy 
c 
b 
2.) 
O-us 
— OZ dy. (10) 
. C2” 
s is identi with equation (9), as may be shown by substituting for 
equation (1 For 50 to be zero it is sufficient for (67u/¢ z*) 
| to be ze odd function of y. The latter does not satisfy the physical 
| litions of symmetry about the y-axis implied in equation (2). Examina- 
{ f equations (1) and (2) leads to the result that the condition 6° 0 
tisfied by solution of the form 
| dU’ 
(y- b*)e i (11) 
dz 
ded that dU, /dz is either linear in z or a constant. While this solution 
tisfies the boundary conditions y% = 0 at y -b, it does not satisfy the 
| itions that uv 0 at two values of z. It may be concluded, 
erelore, that the conditions for dl’ 0 are too restrictive to be met in 
practi 1S 
In Part I th particular assumption for the blade forces led to the result 
| I . . . . a 
7 0 when a, ’ \ similar conclusion was reached in (5). 
tlowever, this conclusion is not of much practical interest on account of 
specialized character of the blade forces. 


5. Trailing circulation outside a wall boundary layer 
In the experiments described, the non-uniformity in the approaching 
w extended to the midspan position. More commonly the flow is only 
uniforn 1 boundary layer on the walls. In such a case vortex 
l 
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filaments giving rise to distributed secondary vorticity and trailing fila.{ 4 TSU 
ment vorticity appear only within the boundary layer region. Owing to the} assum 
distributed secondary vorticity in the blade passage, secondary vorticities The 
are, however, induced outside the boundary layer. Across the boundaries along 
to the ‘passages’ formed by the sheets of fluid leaving the trailing edges fusela. 
a change in w (Aw) will oceur, and these trailing sheets will, therefore | wall ¢ 
contain vorticity. Since there is no trailing filament vorticity outside the | and ¢ 
wall boundary layer, the only component of trailing circulation present | 29. 
must be shed circulation due to a change in circulation about the blades, | gener 

Experiments with the same turbine cascade and a thin approaching | layer 
boundary layer on the walls have shown that the measured circulatior | Inj 
varies outside the wall boundary layer. Good agreement was obtained hence 
between measurements of change in the blade circulation and experimenta layer 
values of Aw in this region of uniform approaching velocity. It was noted | Vorte: 
that in this region, as the boundary layer was approached, C,, decreased with | 
indicating the existence of an apparent ‘downwash’ outside the boundary If. 


layer. As the boundary layer was entered, C,, increased and the major z ¥! 


effects noted in the previous section were reproduced. lead 
upW a 
6. Isolated aerofoils absel 


Isolated aerofoils may be treated by the same method as cascades by| 7. C 
calculating, as shown in (4), the distributed secondary vorticity and, As 
thence, the secondary flows in regions on either side of the sheet of fluid | isolat 
leaving the trailing edge. Values of Aw and the trailing filament circula- | Conv! 
tion (equation (32), Part I) can then be estimated and the shed circulation | filam 
obtained by subtraction as before. Difficulties may arise on aerofoils with | traili 





well-rounded leading edges owing to the inadequacy of theoretical methods | th 
for calculating the distributed secondary vorticity (4). traili 
The strength of the trailing vortex sheet can be estimated approximately | turbi 


from equation (11) of Part I. This is | stren 
the 1 

— A | 2(V/ grad (po/p)).(V .grad)V ; ih aiv(> = (12)} @ con 

q q' q° 7 for t 

TI 


The last two terms on the right-hand side represent the trailing filament 


and trailing shed circulation respectively. If the aerofoil is regarded a aye 
a lifting line with lift Z and approaching velocity U varying in a spanwis tion 
direction, the value of the last two terms, giving Aw, the strength of th TI 
trailing vortex sheet, is the | 


LipdU d & : 1 dL a stall 
~ U2 dz dz cal , 


Aw = pU “pU - : flow 
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result whit h is different from that obtained in Part I, in which JL is 
ssumed t vary as l .. 
The behaviour of the flow about an aerofoil at its junction with a wall 


ong which a boundary layer exists, e.g. at the junction of a wing and 


fuselage, can be discussed qualitatively in terms of equation (13). At the 


| consideration of the secondary flow shows that Aw will tend to zero 


nd consequently the slope of the curve of lift against span will tend to be 


ro. If U goes to zero, then d?Z/dz? will also be zero. Since Aw will 
generally be of negative sign (4), L will tend to decrease as the boundary 
uyer ind wall are approached. 

In general, for thin wall boundary layers, the pressure distribution, and 


Us wil 


ilament 


rded 


panwis 








ence the force on the aerofoil, can only vary slightly in the boundary 
wer region ; consequently, from equation (13), the strength of the trailing 
vortex sheet, Aw, will be small. This conclusion is in general agreement 
with that reached in (4) 

If L is relatively constant and U decreases in the boundary layer, then 

will increase and there will be the equivalent of an upwash which may 
id to stalling as in the turbine cascade. For the isolated aerofoil the 
pwash will be accompanied by a negative induced drag, at least in the 
ibsence of stalling 
7, Conclusions 

\s a vortex filament in a non-uniform flow approaching a cascade, or 
solated aerofoil, passes the aerofoil, the particles of the filament on the 
onvex surface move ahead of those on the concave surface so that the 
filament has to stretch in the wake. This gives rise to a component of 
trailing vorticity, the trailing filament vorticity, which should be added 
to the component due to the change of circulation along the blades, the 
trailing shed circulation. The analysis of some experimental results on a 
turbine cascade with a non-uniform approaching flow shows that the 
strength of the trailing vortex sheet cannot be accounted for solely by 
he measured change in circulation about the blade. The postulation of 
component due to the trailing filament vorticity satisfactorily accounts 
lor the disere pancy. 

The change in circulation along the blade can be estimated satisfactorily 
irom calculations of the secondary flows based on a ‘passage’ theory and 
the determination of the trailing filament vorticity from a simple perturba- 


tion theory 


the results show that the secondary flow influences the behaviour of 
the blade in such a way as to tend to cause a stall near the wall. Similar 
stalling tendencies are predicted for an isolated aerofoil in a non-uniform 


flow and the existence of a negative induced drag is suggested. 
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ON THE DETERMINATION OF THE FLUTTER FORCES 
ON WINGS WITH SUPERSONIC LEADING EDGES 
By B. A. HUNN (1204, Richmond Hill, Surrey) 

d 18 February 1954. Revise received 27 July 1954] 


SUMMARY 


rned pl rily with the practical problem of determining the 


Ly on thin wings vibrating in a supersonic main stream. The class 
f swept wings considered is h that both the leading and trailing edges are straight 
su ; nd the ti traight in the line of flight. Two methods are con- 
to tl iting labour necessary to calculate actual flutter 
s. 7 first is the « | ensive treatment valid for all frequency parameters 

| r small values of the frequency. 


1. Introduction 
é iation of the aerodynamic forces on lifting surfaces due to their 
des of vibration isa problem of considerable importance in relation to the 
ritical flutter speeds. The classical method of performing such 
calculation is to regard the lifting surfaces as of such an aspect ratio that 
two-dime malaerodynamic theory is adequate. This attitude is perhaps 
erstandable in view of the complexity of the integral equations which 
derived from a linearized theory, but it is hardly to be justified for the 


ispect ratio wings now in use on modern high-speed aircraft. 


| lhe problem considered here is the evaluation of these forces on wings 


spect ratio which are swept, and have straight leading and trailing 

iges and tips straight in the line of flight. The basic theory is the linearized 
lynami e valid for wings with supersonic leading edges, although 

d be equally valid for control surfaces with supersonic hinge lines. 
ethods of approach are considered, one valid for all frequency 
meters, the other « for small values. The former employs a trans- 
iables w] permits the ultimate numerical evaluation of 
nd the latter employs a set of potential functions derived 
vance and t vuthor (unpublished work) for a prescribed 


+ 


the class of wings mentioned. 


aetormation mode $s oO 
It has been shown by ¢ ck and Rubinow (1) that the potential ¢,- of 


ipper surf fa wing which is in a truly supersonic stream 


i from 


L [| eet) &(t—72)} dody (1) 
2a J. J{(~a—£)?— B2(y— n)*} 


f ed aerodynamic theory. In this integral the point (x, y, 0) 
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is the one under observation and the point (€, 7) is a current one which may} © tt 
traverse the area of the wing in the forward Mach cone from (2, 7,0). Thy the 
upwash velocity w(x, y, t) on the wing is assumed to be a separable functio, 


of the form wha) W(x, yet) 





It is: 
and the intervals 7, and 7, are given respectively by | force 
M (g—é) 4&2 £)?— B*(y—7)? It is 
T 2 “79 A S J : 9 (3 a6 
Bea, Ba, the ¢ 
with B? M*—1, where M is the free stream Mach number and a, is th 
free stream speed of sound. 
The local perturbation pressure on the upper surface is given by 2.7 
W 
| l Qby , by 
Pr pV — ——-+- — { vote 
V ot Cx 
and the local lift per unit area by Oy 
) 
I PL PU» 
p,, being the pressure on the lower surface. 
Finally, the section lift and moment for a wing are given by lian 
T.E. 
lift/unit leneth | (Pr Pr ) da, { 
L.E. 
rE Bi 
moment about leading edge/unit length | (w@—ay,..)(pp—Ppr) a. sess 
L.E. “7 are) 
Since, however, within the framework of thin plate theory, the perturbatio 
pressures on the upper and lower surfaces are of equal magnitude but} whe 
opposite sign, it is sufficient to consider equation (6) in the form aspe 
T.E. 
201 l éd,, . Ody ii (8 
2p = a +) ae, 
J VW at Ca 
L.E. the 
with a similar result for the moment about the leading edge. 
When the motion of the wing in the stream is simple harmonic, It ! db, 
convenient to write -" 
@(t) pial 
and dr, (X, ¥, 0, t) ,-(x, yeh. (I 
It is 


It then follows that 


- iu od - 
ol ca eb +e } ( 


OX 
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that the determination of the section lift and moment depends upon 
the quantities l.E T.E. 
rE Y | dy (x, y) dx; |} rdp(x,y) dx. (12) 
‘T mee 


[t is evident therefore that the determination of the unsteady aerodynamic 


es on a vibrating wing can be made to depend upon the flow potential. 
s sin y true for modes of vibration having chord-wise distortion, for 
dditional strip derivatives will depend upon 


UPL Pr ) dx. 


2. The evaluation of the integral for the potential 








When the motion of the wing is oscillatory, the integral defining the 
haved be written in the form 
[ [ W(€, 7) 
res « 2_ By 7)?! 
w,/1(x—&)?— B?(y— n)?! 2s : 
OS} \ —____— ! > ad iK(a $ dédy (13) 
3A 
e eg 10 where A is defined by 
Mw 
K —_ (14) 
Ba, 
Before dering the means of evaluating this integral, it will be 
venient toremove the scale effect. For this, non-dimensional coordinates 
ntl 1. defined by . - 
] ined i X alé. 
) y/b = y/Aée, (15) 
geometric mean chord, 6 the total wing span, and A the 
spect rat Modifying the definition for (€, 7) and using 
S? (X —&)?— B?A*(Y — n)?, (16) 
nteg ecome 
l i ._M? oo. a S 
I pI i “of (X—£) — 
S B M 
i ( + psy Sew 
exp 1——y * Gat dédy. (17) 
P| B Ta) | 
is convenient to split this integral into its two parts: 


6 (ry .M? ; Bis .. 
exp! ‘= (A é)- } dédy (18) 
at) Ss | B 


M 
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. brrW_i .Mt@/~ . S\\,, 
and db}? — 5 | | zr exp| —f Re (a —£)4+ a déd». (19 
Then, using the transformation 
ee ‘ " ‘ : 
w= Spvi(M +1) MiL V(X —€)+ BA(Y —n)}—y{(X—£)—BA(Y—»)}], 
te 
v= sR. ((M —1)vMj[,/{(X —é)+ BA(Y — n)}+,/{(X —£)— BA(} n)}l, 
(20 
dy) becomes ' bP = | | We-t+e) dudv, (21 
and using the transformation } 
l : . " 
‘= api 1) Mv}| V(X —€)+ BA(Y — n)}+,/{(X —é)— BA(Y—7) 
v sp vie —1)Mrj[y{(X €)+ BA(Y —n)}—.,/{(X —é)— BA(Y —7)}], 


»)») 
4 


the expression for 4%’ becomes identical in form with that given by (21), 

It should be noted particularly that points within the forward Mach 
cone transform from (x, y) on the wing, to points in the (wu, v)-plane withir 
the sector bounded by the transformed characteristics in the first an 
second quadrants in the first case, and in the first and fourth quadrants 
in the second case. Finally, writing 


rcos@ u| 
rsin 6 » J 


both integrals reduce to the form 





bu = =u || W(r, 0)r e~*”* drdd. 24 


3. Transformation of the boundaries 


(i) The transformation of the characteristics through (X,Y) given by 


(X—£) = + BA(Y—7) (25 
can be shown to give, for each transformation, the equation 
= en Ja \° = 0 26 
\V(M+1)* (a1) 


These are repeated lines whose inclinations to the u-axis are given by 


M—1 


} F ae | . 
alias F/G) 
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(ii) Characteristics through the point (Xy, ¥,) may similarly be seen to 
transform into 

j u 12 ' F ; , ; : My 

' (Xx —IJ+ RAF —T)i- 

\,/(M 1) /(M 1)} I 0) 0) 53 


For each sign there are two lines, one inside the appropriate sector of the 


(28) 


)-plane and one outside. This fact will be needed when tip corrections 
wre to be considered. The point (XQ, Y,) lies on the given characteristic. 


iii) The transformation of the wing leading edges may be written in 


the form mé +An. (29) 


Again, manipulation of the transformations shows in both cases that the 


leading edge become 


») ‘ 
mus wie thi ais , Mv ‘ 
4 (mX -AY)—.. (30) 
M+1° M—1— B B 
Since ot A, A being the leading edge sweep-back angle, and since 
furthe1 , ? 
Bin l (31) 
for supersonic leading edges, equation (30) will be identified as a pair of 
ellipses centred at the origin of the (w,v)-plane with major axes inclined to 
the LXI1S noles 


I ? 
btan ies ). (32) 
m 


The relatic between ind 6 can be seen to be 


vM sin A(mX FAY) 


Sal sane Te Oe (33) 
V/ COS \ cos(20-+-A) 
The upper sign is appropriate to the starboard leading edge, the lower to 
the port 
Finally, using the plane polar coordinates, the @-coordinates of the 
ntersections of the two ellipses are given by 
AYM 
! 2 weer (34) 
,/(X?+ A?Y?) 
1 y ~ 
vhere tan-* — tan-!: (395) 
- t 
4. The transformed upwash function 
For flutter purposes it convenient to regard the wing as a lamina. 
A | 4 } 
ANV distor ! may be represented by 
, — = Ap, x" y?. (36) 


f the distortion shapes, or, for that matter, 
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their experimental determination, justifies limiting this series to a fey 
terms. 

The upwash condition can be shown to be 


C = tf 
— (x, y,0,t) = w(x, y, t) V—+—, (37 
Cz ‘ ; Ox «ol 


where z g(x, y,t) describes the wing motion. For the oscillatory case, 


separable functions can be used to give 


c am G a y 
— ole, 9.0) J LiwG, (38 
C2 : Ox 

where g(x,y, t) = e&“G(a, y). (39 


Observing that the series form is still valid, one finds that 


— 4 2 Dp 2 7 
S33 a $? 7" SN “a x am (.W—cos 29) ly ; 7 sin 24 , (40 
a pq > hat Lag PH U ] I Al 


where a,,, is in general complex. 
Using this series for the surface configuration of the wing, it then follows 


immediately that W — Wor.6). 4] 


2 


In this event, since | pen ye-it* dr 


is immediately integrable, it is convenient to perform the integration with 


respect to 7 before @ in the determination of ¢,,. 


5. Zones of integration 
It can be seen without difficulty that the total potential on the upper 


surface of the wing can be represented as 


0; RO) 
IG a a —— oo 
oy (x,y) = dP +42 = — iT; dé | W(r?, O)re-”™ dr. (42 
TTLLV 
6" 0 
where 
; /{M—1 
0 tan-! l(a 
Vv * (43 
[ (i ] 
Ge aw—tan-! / . 
J \M-+1, 
and 7 R(6) is the sole means of discriminating between whether th 


point (2, y) is forward or aft of the Mach cone from the leading edge vertex 
and whether this potential is subject to tip effects. 

For points on the starboard wing forward of the Mach cone and 
not subject to tip effects, the zone of integration in the (w,~)-plane 
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I responding to the interior of the forward Mach cone from (a, y) occupied 
by the wing is, in fact, that segment of the ellipse 
nu p- 2uv Mv , , 
: == ——(mX —AY) (44) 
. M -+-1 uf —] Bb Bb 
if 7 7” 
ine —_— 0) (45) 
M+1)' /(M—1) 
0 first, second, and fourth quadrants (see Fig. 1). 
‘or point ing within the Mach cone, again not subject to tip effects, 
e regi f integration is bounded by two ellipses. If, on the starboard 
le. t} teo} s taken over the range 
{Aj | 
ta | . ) a yb, 
J \M-+1/) 
a * (46) 
. pM sin A(mX+AY) 
M cos A cos(2¢ \) 
0 us, 
vM sin A(mX —AY - 
' LA ) (47) 
| MeosA cos(26 \) 
{{M —1 
7—tan-! | } (48) 
Aj \M-+-1 
g inequality (46) and where y%, and %, are the two roots 
34) and (35) lying in the half-plane to the right of the line given by 
see | 
Fj ints subject to tip effects, use is made of the results given 
VE rd (2). In the (a plane the area of integration is limited by the 
| lorward Mach cone from /), the characteristic from that point on the 
| edge at the tip intersected by the complementary characteristic from (x, y) 
Ul ( ng edge 
he esponding zone in the (u,7)-plane is either of the two cases 
previously mentioned but with a segment sliced off by the line given by 
: )? : . , ... My 
| _R Q) 
—— 7 ait [((X—X,)—BA(Y—Y,)] RP" (49) 
In parti t should be observed that if Y, is evaluated at the tip and 
nding va of € for the intersection of the line 
\ BA(Y¥ —yn) (50) 
Y= 9, (51) 
roaches the tip, the nearer the line given by 
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(49) approaches that given by (45). Finally, when coincidence occurs the 
region of integration is zero (see Fig. 3). 

This treatment of the tip correction needs some modification in the light 
of information received since the original draft of the paper was prepared, 
It had been stated by Stewart and Li (3) that the Evvard method would 
be suitable for oscillatory motion. This has subsequently been retracted (4 
and they point out that errors of the second order in frequency would be 
introduced. Another method has been given by Stewartson (5). 

No attempt is made in this paper to utilize their suggested method of 
correction. It should, however, be borne in mind that flutter usually occurs 
on main surfaces for frequency parameters well below unity so that th 
errors incurred by neglecting these effects will probably be small and 
restricted to small regions, but control surface flutter, occurring for values 
of frequency parameter usually greater than unity, will be affected more 


seriously. 


6. The low-frequency approach 
In this case the aerodynamic theory is again linearized with respect to 
the frequency parameter, the assumption being made that only first powers 
are significant. On writing 
M2 x 
p(x, y, 2, t) D(x, y, z)exp iw| t— — - 52 
QD\4, 4, 6, HY, «Je, WwW ke - 
a i 
with the previous hypothesis, it can be shown that the integral defining 
the potential becomes 
brrw . M? 
O(z, y,0) = —— | ay—— & } dé 53 
D(x, y, 0) - | <q La ef) ten (53 
using the non-dimensional variables. 
An equivalent approach from the equations of fluid flow shows that ® 
must be a solution of 


2) 2D 2) 
— + 9 (M?- . os if (94 
Cy2 az ox" 
Using characteristic coordinates of the form 
r = «— By) (es 
s=a+By)’ 
Ward has shown (6) that ® is given by 
.. Fe Ww’ , 
® -—— | | saieroal imatensen dr, dsy. (ot 
27B JJ vi(r—r.)(s—s.)} 
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304 B. A. HUNN 
Here W’ has a slightly different connotation, being in fact derived from 
e® 
~€ 

C2 
or, consistent with the low-frequency parameter approach, 


9 \ 


alt) 


. M*2x 
W’ = Wie, (1 Wn | 58 
be (v7, y) w RB } ( 


\ 
It is not difficult to show that (56) and (53) are identical. The form of 
(56) is, however, much more convenient for the determination of analytical 
expressions for the potential. This form has been used to compute th: 


potentials for wings whose incidence is given by 


C 
oe xa y* 
Cx . 
(59) 
C2 
and vw" |y® sgn(y) 
Ox 


for all positive and zero values of 7 and s such that 


. — » 
'—so ”. 


These include tip corrections as required for the class of wings defined in 
the summary. The analytical expressions obtained have been written in 


terms of the two parameters 


a Bm (see equation (31)) 
d 
and 7) BY. (60 
xv 


Some functions have been computed for a few values of a over the corre- 
sponding ranges of 7. Those omitted have necessarily been the ones con- 
cerned with tip effects, for they depend upon the geometry of the specifi‘ 
wing. 

Using the reduced potential it is convenient to work in terms of the 


components defined by 


@ 0,+09, (61 
where 
aty) em ,OG \ 
1— W, J 
Cz Cx 
(62 
2 IG M22x0aG 
oo, = W;, Ae Tpsoalle 
a a Cc B- C £2 


with G defined as in (39). 
Further, the section lift can be shown to depend upon 


T.E 


(2 j i. ues) dy (63 


te C Ox 


\ca M? 


iv(M?/B2yr/c W (x,y), (57 
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rom | witha similar expression for the moment. It then follows that the in-phase 
depends simply upon ©®,(xy.,,,¥) while the out-of-phase lift depends 
rE. T.E 
| r 1D) F = o® 
D(X 4. Y =} — dx | M?~— —1 dz | (64) 
fp” J c P C CL 
oy | LJ L.E ; 
} 
Finally, it should be observed that the reduced potential can be compared 
with the potential given by the other method for low values of frequency 
meter O the basis tl 
@ ®,){1—a } ~, inf ,. ~~ YD, b7;>(2, y) (65) 
=a B- si B? Cc ; 

7. Numerical comparison of results 
\ swept wing having 44° of sweep-back on the leading edge was con- 
The modes of vibration prescribed for this wing were uniform 
unging motion and pitch about an axis through the leading edge vertex. 

In the first se, then 
Lu re & » a 
W(x, y) Viv, (66) 
F 
| being the amplitude of oscillation. In the second case 


i / Vo(1 iv). (67) 


For each of these modes the potential at two points on the same section 
the wing was calculated, one outside the vertex Mach cone and the 
. ther insid 
The integ can be shown to reduce to 
each 
| | {ve(cos R2—1) —°°8 =" sin R2— R? cos (1 7): 
7™Myv J | M M /i 


spas sin R?+ R2sin RI ral dd (68) 


M MM 
he pitching case, R(#) being defined by 
VE «i YAY 
R2(6 vd sinA(mX FAY) (69) 
VU COS \ cos(26 + A) 
the upper sign be ing taken for the point outside the Mach cone on the 


xX 
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starboard wing and the alternating signs as defined in equations (46)-(48) 
The limits 6, and @, are defined in (43). For the plunging case the integra] 
becomes 


2 





dy = — Ve - - | {(1—cos R?)+isin R?} dé (70 
C 7ihv ; 
01 } 


with R(@) and the bounds of @ defined as above. 





Consideration of this integral shows that the chances of deriving analyti- 
cal results is small. A numerical integration was carried out using the 


generalized Simpson’s rule. In particular for the point outside the Mad 
cone the integrand was evaluated for nine different values of 6, and for the 
point inside the numbers of values of @ taken were first three in the range 
defined by (46), then eleven in the range defined by (47), and finally three in 
the range defined by (48). The values of v which were adopted were 
vy = 0, 0-1, 0-2, 0-8, and 1-4. 
The other parameters were: 
M 1-84, 
Ar = 137%, 
X 1-9266 inside; 1-3994 outside, 
Cc 10-41 ft. 
The rather peculiar values were chosen so as to permit comparison wit 
a previous calculation using the low-frequency parameter theory. In par 
ticular the value of a for this aircraft at the given Mach number was | 
These results were then plotted against frequency parameter and the 





corresponding low-frequency parameter results were superimposed. These 
had been obtained direct from the graphs mentioned in the previous 
paragraph. The slight disagreement between the two methods at » 

should be observed as due to the inaccuracies in graph reading and thi 


inherent errors in numerical integration (see Figs. 4—7). | 


8. Discussion of results 

The criteria which must be applied to any numerical technique requilt 
for practical work are that it possesses a logical simplicity, that 
involves relatively simple numerical operations, and that the volume 
operations is not prohibitive. With the compatible frequency paramett! 
theory it is evident that not only must the complex potentials be evaluat 
at a fairly large number of points spread over the wing, but the integ! 


involved for each point must be repeated for each new value ot ' 
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fequency parameter. Further, for every new mode of vibration these 
on ust be repeated. It therefore follows that the volume of 
involved in this comprehensive approach is large. 
On the other hand, the low-frequency theory can be used for all reasonably 
lues of the frequency parameter, and the numerical values, except 
re tip corrections are involved, can be obtained from the functions 
lv tabulated. The groundwork which has been done on tabulating 
ene a in be regarded as common to all calculations involving the 
ss of wings previously defined. From the numerical point of view there 
therefore much to commend this method. However, from the mathe- 
Pare tical point of view, the hy pothesis that 
j . Mea . Mx 71 
exp| U1 Bee ~1l—iu RB = (71) 
sufficiently small values of v would seem to be a very restrictive one. 
Now it usually happens that main surface flutter occurs for values of 


juency parameter not greater than 0-5. Referring to the graphs for the 


ntial it will be seen that the discrepancies between the two theories 


large for 0-5 except in Fig. 7. This, however, must be considered 
he main stream Mach number, for the basic integral involves 

frequency parameter in the forms M*y/B and Mv/B, and one would 
rreater discrepancies between the two theories when the wings were 


v just supersonic. Bearing this in mind it can, however, be seen that 


ne proceeds with due caution some useful results can be obtained from 


"WV irequen¢ \ theory 


$$ 
J ~ 


\ complete flutter calculation has been performed on the wing under 
sideration, using the low-frequency parameter theory for comparison 
th results predicted by incompressible theory. Whereas the result pre- 


©. licted for incompressible theory gave a Mach number of 2-8, the supersonic 
20089 | theory gave M = 4-2. This result thus appears to be rather dubious. 
However, it has been shown experimentally that the ratio of the experi- 

ital wing flutter speed to that calculated on a theoretical, incompressible 


isis in the supersonic régime increases considerably above unity with 


oo 


—— 
Mach number (7 


Finally, although specific mention has been made of wing flutter, the 
methods discussed ean be ipplied to control surfaces. The chief disadvan- 
tage occurs using the low-frequency theory, for one usually finds frequency 

rameters greater than unity arising from control surface flutter. There- 
lore it would be desirable to use a third-order theory in frequency as devised 

Nelson (8) as an alternative; but the more comprehensive method 


lescribed in the earlier section of this paper should prove the more accurate. 
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ON THE BOUNDS OF EIGENVALUES 


By K. WASHIZU 


(Department of Aeronautics, University of Tokyo, Tokyo, Japan) 
Received 23 September 1954] 
SUMMARY 
n this formulae for locating the bounds of eigenvalues of a real symmetric 
Vy} A 
itrix are treated in a unified manner. In the first part of the paper, several theorems 
which provide bounds for eigenvalues are represented in a geometrical space. In 
, second part, a generalization of Rayleigh’s principle is attempted, and some 
ems which provide bounds for eigenvalues are discussed as special cases of the 
S) | 
generalized Rayleigh’s princiy The system treated here has a finite number of 
A ugus 8 I i 
grees of freedom. By appropriate modifications the above principles can, how- 
r r, De xtend itoa vstem with an infinite number of degrees of freedom. 
VW 


1. Introduction 
WHEN an eigenvalue problem can be solved only approximately, it is 
lesirable to obtain bounds for the eigenvalues. Several theorems pertain- 
g to such bounds have appeared in the literature (1). Among them 
Rayleigh’s principle for estimating natural frequencies of vibrating systems 
is undoubtedly the most well established. Rayleigh’s theorem shows that 
his approximate method always provides an upper bound for the lowest 
irequency and a lower bound for the highest one, and that first-order 
rrors in the assumed mode of vibration produce only second-order errors 
the estimate of the natural frequency. In practice a criterion of this 


] 


ind is of limited value unless it is accompanied by another theorem 





vhich provides a lower bound for the lowest frequency or an upper bound 


tor the highest one 

A formula giving a lower bound for the lowest eigenvalue was obtained 
by Temple (2, 3). Kohn (4) and Kato (5) generalized his theorem so as 
to determine the location of any eigenvalue, i.e. not necessarily the lowest 
tthe highest one. Later, Temple (6) showed that the Kohn—Kato theorem 
n be extended to a system with inertial couplings, while the author (7) 
owed that these theorems have definite geometrical meanings and 
represented them in a vector space. A brief summary of the principle of 

geometrical representation is given in the first part of this paper. 


Recently R. V. Southwell (8) published a paper giving new theorems 
-™ , > » ° . 
the location of the bounds of eigenvalues. His formulae seem to 

suggest a method of generalizing Rayleigh’s principle and have stimula- 


ed the author to write the present paper. In the second part of this 


{Quart. Journ. Mech, and Applied Math., Vol. VIII, Pt. 3 (1955)] 
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paper a generalization of Rayleigh’s principle is attempted and the theo. 


rems mentioned above are discussed as special cases of the generalized 


Rayleigh’s principle so that their mutual relations might be understood 
clearly. New formulae could be obtained from the generalized Rayleigh’s 


principle, but not all, probably, will be of value for numerical computa. | 


tion. 


2. Rayleigh’s principle 
Let us consider an eigenvalue problem of the type 
Au Au, (] 
where A is a real symmetric matrix and wu is an n-dimensional vector wit} 


components 
ponents (084, ay-<5 Hq) (2 


Let the normal modes or eigenvectors of A be denoted by 4,, 4,,..., 
and the corresponding eigenvalues by Aj, As,..., A,,, so that 
Ab, = A464, OF = 1,2,...,%). 3 


All of these eigenvalues are real. Therefore, let the order of magnitud 
of the eigenvalues be specified so that 
hy Me Ku. HAY 4 


The normal modes are orthogonal and normalized in the sense that 


(¢,, db, ) Onc 5 


rs? 
where 6,, is Kronecker’s delta, i.e. 
Sng -0 ifrss, 
5 l tftr=e. 


We shall henceforth choose the » orthogonal and normalized vectors as 
the bases of the coordinates in the n-dimensional vector space and let ' 


be an arbitrary vector which has components of (¢,, ¢s,...,¢,,) with respect 
to the base vectors ¢,, so that 
n 
wo = 2 or dr 6 


With the above definition, the following theorem, referred to as Ray- 
leigh’s principle in this paper, is readily obtained: 


tayleigh’s principle. 

< (wy Aw) — 
(w, w) 

The proof of this theorem is easy and need not be given here ; the theorer 

holds irrespective of the signs of the eigenvalues. 
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theo. For the prol lem of the type 
ulize Vv = AT 2, (8) 
namely, a system having both elastic and inertial coupling, the theorem 
c ; ‘ 1 . . cal . . 
expressed in (7) can be extended with a slight modification. It is assumed 
t V is a real symmetric matrix and that 7' is a positive-definite real 
} ma ? 
symmetric matrix. Now, by putting 
‘es he by A, Tv u, (9) 
8) can be transformed to a type which is equivalent to (1). Therefore the 
responding relations can be obtained immediately. Thus, in place of 


Vey, A, Tys,., (10) 
T's. ) O, . 
| | (11) 
Vb.) A,.5,¢ J 
nn 
ue SN co, (12) 
r—1 
| | 9) ee t 7) we et 
(w, Vw 
ee ek (13) 


3. Vaisey’s theorem 
Southwell (8) states Vaisey’s theorem in the following manner. For a 


tinuous system, consider a differential equation of the form 


i 
i boys vow 0. (14) 
ors a | where A of Southwell’s paper is here replaced by v. Let # and # be self- 
1djoint linear operators and let 
. | , 
| l | ww dudy, (15) 
I /A\ || wi'w dady, (16) 


Ray- | lere it is assumed that / is positive definite and I’/A may have any 
sign. As a consequence, the eigenvalues v can have both signs, but they 
must still be real. Let the positive eigenvalues be denoted by 
’ 19 Var-ees Vi 
nd the negative ones »\ 
1 9 V l 
with both series in order of ascending numerical magnitude. With the 


initions, Vaisey’s theorem may be stated as follows: 
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VAISEY’S THEOREM. The value of vp defined by 


AJ | | witw dady 
VR Ty’ a : ( 17) 
| | ww dady 
cannot lie between the least positive and the largest negative eigenvalue (; 
and v_,). 

Now it is easy to show that there is a close similarity between (8) and 
(14) by writing V, 7’, and A in place of 3’, 3, 1/v, respectively. Thus, the 
conditions which define each problem become equivalent, and Vaisey’s 
theorem leads to the statement that the value of A; defined by 

(w, Vw) 
Ar —, (18 
(w, T'w) 
can only lie between the least and the largest eigenvalues (A, and A,) if 
we remember that 1/y, and 1/v_, correspond to A,, and A,, respectively 
Thus we see that Vaisey’s theorem is contained in (13). 

For the sake of simplicity, we shall henceforth treat only problems of 

the type (1) unless special reference is made to problems of the type (8). 


4. Principle of geometrical representation 

We shall consider an eigenvalue problem of the type (1), and let w by 
an arbitrary, normalized vector as expressed by (6). Now, consider tw 
vectors w (OA) and Aw (OB) in the n-dimensional vector space as show! 
in Fig. 1. Points A, &....; and P, are located on the line OA so that 














OP, = A,. (19 
B 
7 
wy / 
0 all 3 NV 
w PCB P., 2B 
/ 
Fria. 1 


Since 








BP, . BP, cos ZP, BP, = (BP,, BP,,,) = (A;w—Aw,d,,, w—Aw) 


i 


: >) (A;—A,)(A; +] —A,)e 


r=1 


it follows that LP, BP; < 4. (21 
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| 








Now, 
OAB 


is the 


or eq 


Fron 
Le 
let th 


the n 


case 


and 
Thes 
who: 
A, re 
the 
be e: 
be 1 
relat 
5, k 

O 

L 


and 


T he 


bh is 








ON THE BOUNDS OF EIGENVALUES 315 


Now, let BC be a perpendicular from the point B to the line OA in the 
0A B-plane Then 

















17) | a 4 
7 (Aw, w) (OB, OA) OC (22) 
is the Rayleigh quotient. By using relation (7) it follows that 
P, lies in (—, 7) 
l | se 
sage ; (23) 
| and P, lies in (7, +9) } 
ss or equivalently ( 10) BP, < hr ) 
; (24) 
ISe} und 2 P, B(+ 00) : la) 
From (21) and (24) the following conclusion is derived: 
Is Let M and N be two points on an arbitrary, normalized vector w OA and 
\ the point B be defined by the re lation Aw OB (Fig. i. If 
MBN > }r, (25) 
at least one eigenvalue must be contained in the interval MN, i.e. 
mS OL | OM <A, < ON. (26) 
\ very simple example of this relation is illustrated in Fig. 2. In the 
ise 7 yA we obtain 
cos 0, sind >< 
A 0\ /eos 6 A, cos 8 \ 
Pend, el, 
7” 0 A,/\siné@ A, sin @ B 
P, BP. lz (from (20)). } 
' : , / 
These relations are used in drawing an ellipse / 
whose major and minor semi-axes are A, and . r, 
respectively To the problem of type (8) —— 
° a" ‘ 2 
the above geometrical representation (26) can Pe 


e extended with a slight modification. By using (9), the problem (8) can 
reduced to the same type as (1), and, in this modified space, the 


ution (26) holds similarly 


5. Kato’s theorem 


One of Kato’s theorems (5) can be stated as follows: 


Let be an arbitrary ywmalized vector, so that (w,w) 1, and let ”) 
] € h, li I as follou ‘ 
fa. w) (Aw— nw, Aw— yw) e*. 
Then. for « r such that a , (a,b) contains a point of the spectrum, where 
ij ] 


27) 
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Now we find in Fig. 3 that 

n = (Au,w) (OB, OA) OC 

e? = (Aw— nw, Aw— nw) BC? -- (28 


and e*+ 7? (Aw, Aw) 


line OA such that a = OF < OC and 
EBF = }3n. Then we get 
EC.CF BC? = e, 
\ rp and 
OF — OCL0F = 94-2 ik 


1) a 











This is the geometrical meaning of 





(27), which is already evident 
from relation (26) if the conditior 
EBF = 37m is taken into consideration. 
6. Temple’s theorem 
Temple (6) extended Kato’s theorem to the problem with the type 
Vv = AT», (8 
namely, to a system with inertial couplings. His theorem may be stated 
as follows: 


Let w be an arbitrary normalized vector, so that (w, Tw) 1, and let p an 
€ be defined such that p = (w, Vw) 
R= Vw pT'w 3 (29 


e} = (R, TR) | 
Then, for every a such that a < p, (a,b) contains a point of the spectrun 
where b is defined by (b- 


. ‘ 
-p)(p—a) €p. (30) 
Temple’s theorem can be represented as in Fig. 4. The geometrical 
RB 




















relationships are similar to those for Kato’s theorem. 
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7. Kohn’s theorem 
The main parts of Kohn’s theorem (4) can be stated as follows: 


P Some eigenvalue A, lies in the interval, 
Fi {2 9 2|% « < A2 9 2) 
1*)av— 20Ag, . Aj oa [(A av 20 Ay - It. 
(31) 
} here S wbitrary ve number, 
A.. Aw) and (A?),, (w, Aw), 
ind a wbitrary norma ized vector. 
Now. by transforming (31) into the present notation, we get 
|e7-+-(a 2]2 <A; < at+fe?+(a—7y)? I]t. (32) 


Let us consider the geometrical meaning of (32). Take a point G on the 
e OA so that OG ., as Shown in Fig. 5, and let H and F be two points 
the same line on opposite sides of G so that 

: EG BG GF. (33) 

Then the relation (32) states that 


OE <A; < OF. 

















he equi ence of Figs. 3 and 5 is then easily seen. 
Be 
— 7), ———>} | 
| —— a — 
; | Fic. 5 
| \ GENERALIZATION OF RAYLEIGH’S PRINCIPLE 
| 8. Matrix function 
Consider a rational function in terms of A such that 


( A 
f(r) — 2A) (34) 
h(A) 
x, A P* sea “T Xm A”, 
h(X) = Bo-+By A+ .---E BX, 


are arbitrary real constants. We can build 
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a matrix function to correspond to (34) in terms of a real symmetrical 


matrix A such that 

f(A) = g(A)[h(A)}-? = [h(A)}-19(A), (35 
g(A) = ag I+a,A+...+a,,A™, 
h(A) = ByI+B,A+...+B,A', 


and J is a unit matrix. Then the following theorem is obtained: 


where 


If we denote the eigenvalues of A and f(A) by A, (i 1, 2,...,2), and 


K; (¢ = 1, 2,...,m), respectively, then 


K; f(A;) (7 ee nee (36 


Thus there are correspondences between A; and «x,, although the orders 


of magnitude and degeneracies are in general completely different. 


9. A generalization of Rayleigh’s principle 


Let us now consider how Rayleigh’s principle can be generalized in th 
case when the matrix A is transformed into f(A). For convenience, le! 


the modified Rayleigh quotient be defined by 


a 


d (w,f(A)w) 
nani (w, w) 

Now let the relation « = f(A) be plotted as shown in Fig. 6, and tet thi 
minimum and maximum eigenvalues 0 
[(A) be denoted by «,, and «,, respectively; 
A generalization of Rayleigh’s principle is 
obtained by applying it to the matrix 
f(A), giving 


Kp : Arw : Ky: {od 





Let us denote the abscissae of the points of 
intersection of the line « = Ap.) with th 
curve x = f(A) by a and b.t| Then a geo- 
metrical interpretation of the above rela- 


tion is as follows: 





If Ap, 1s defined by (37) and two points 
and b are defined as above, then at least one eigenvalue is contained intl 
inside and outside of the interval (a,b) respectively: 

: =) 

a A, ) | 

and r, lies in either (—oo,a) or (b, +c) 


For simplicity, we consider a case in which « = Apy and x f(A) have only two int 


section points. It is easy to extend the result (39) to cases with more intersection pol 
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0. Calculation of A,,, 
It is not usually easy to calculate the matrix [h(A)]- in (35). For 
alculating A,,, alone, however, the following transformation may be 


seful. First, let us assume that 
R(A,) (¢ = 1, 2,..., n) (40) 
re all positive. Since we can find a real vector wu such that 
h(A)|-w = u, 
, can be transformed into 


h(A)| lg(A) (w, [h(A)| tq(A)[h(A)]-*w) 7 (u, g(A)u) 


(w. w) (w, w) ~~ (uw, h(A)u) 





writing w in place of u, we obtain 


RM = Gp, h( Ayo)’ 


\ (41) 


Therefore, if assumption (40) holds, then A,,, can be calculated by (41). 
| RM A 
If assumption (40) does not hold, then we transform f(A) into 


h(A)g(A) ) 
f(A) = ——~— > (42) 
| h(a) |? 
nd by applying (41) we obtain 
(h(A )w, g(A)w) ° 
a lc th at Se (43) 
(h(A)w, h(A)w) 
Che two relations (41) and (43) may prove convenient in practical calcula- 
In the following sections it will be shown how the theorems which 
provide the bounds of eigenvalues can be obtained as special cases of this 
generalized R vleigh S prin iple. 
11. Case 1 
vhere in arbitrary real constant. In this case we have 
1) (A—al)-! (44) 





(45) 
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then the eigenvalues of f(A) are arranged in order of their algebra; 
magnitude, namely 





] l l l l 
a — 5] eres 7 ' Teens | eee, —— ~* 
A;— X A; | xX A, x A, on (7 A; 1 X 
Therefore the generalized Rayleigh’s principle leads to 
I l 
eee A S —- —— (4f 
RM , 
A;—« Nj43—0 
(w, (A—al)-lw 
where Any = —— Bi (47 
? (w, w) 


or Arv = = ‘ (48 








>» 
=> 
xs 
> 
s 
a 
—— 
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12. Case 2 
f(A) = (A—a)?, where « is an arbitrary real constant. In this case we 
have f(A) = (A—al)? (49 
and “,= (A,—a}*® (¢ = 1,2....,2). (50 


The generalized Rayleigh’s principle asserts that if A,, and A, are defined 


so thé 2 : 2 
so that Ky = (A,—«)? = min(A;—«)? ) 


—a)? max(A,—a)? } 


~] 


then (A,,—a)? < Ary < (A, 





where Ary = 


These 
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ne These relations are shown graphically in Fig. 8, where « = f(A) and 
\»,, have only two points of intersection, the abscissae of which are 
noted by a and b, respectively. Then we get 
a A hh. (54) yr x) 
} \ If ! eithe ( f t) Ol (bh. 0 }. (55) 
In orde ilculate a 1b, let us define e and 
that 
so th ™ 
* lw, Aw), 
he rmalized condition that (w,w) l. 
Fic. 8 
| we LI? 
e- ”" 2an x" 
€~ {x n)* 
x Vie" (a n)*} Ti 
. (06) 
x Vien (x n)*\ 
s easily seen that relation (54) combined with (56) is equivalent to 
Kohn’s theorem, which was mentioned in section 7. 
13. Case 3 
} Y 2S 1) bitrary real constant. At the outset, let us 
ssume that the eigenvalues of A, namely, A; (¢ == 1, 2..... n), and « are 
tive Is Case We h € 
{ x 14+aA-! (57) 
. | 
A x ° > ~ 
(2 ep ieee nN). (58) 
A, 
All the «; are positive and equal to or greater than 2. The generalized 


hayleigh’s principle asserts that if A,, and A, are defined so that 


(Ac, =| 
\ x A;) 


min 


\ x (A; 


x 
K max — 


| 
A, la § A,J 


: 
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(w, (a-1A +aA~-!)w) 














where Apy = (6] 
_ (w, w) , 
(Aw, Aw)+a?(w, w) 
or Arw = . (62 
x(w, Aw) 
These relations are shown graphically in Fig. 9, and we get 
y < r,, x < b, (63 
and r, v lies in either (0,a) or (b, +), (64 
where . 521 98 
a 1 +-d—.,/(0°-+ 20), 
b 1+-8-+-,/(6?+ 26), (65 
and Ary = 2(14+8) 6>0 
A,@ 
_-_ oe 
Arm 
| ace aces 1 
oe. 
! 
i$ th 
-4 1! ; «8 
; a |i 4 
| 
| 
! 
ee 
Fig. 9 
Next, let us remove the restriction that A, (i 1, 2,...,2) and « area 


positive. Unfortunately no useful criteria for the bounds of the eigen- 
values are obtained in this case, as is easily shown in Fig. 9. For exampl 


even if 5 > 0 still holds, A,,/« can lie either in (—o, 0) or in (a,b). 


14. Case 4 


» A : x 1 xA ° ‘ 1 
f(A) = _ = = whe renwan arbitrary re al constant. In this 
OX A A? -o : 

case the transformed matrix is 
f(A) (a 14 r 1A 1)-1, (66 

and the eigenvalues of f(A) are 
A; a\-1 A; : " 
Ke = (44 ea es. eee (67 

X A; A? oe 


Clearly, f(A) is an odd function and bounded. Therefore we can expect! 


have better formulae for bounds than those given in section 13. Now th 
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seneralized Rayleigh’s principle asserts that if A,, and A, are defined so 














that \ 
X/ . 
K,, —P mint —— | 
‘ A‘ ia \A? 1 x2} . 
(68) 
XA, | x i 
K,, et max} — - 
; AG x“ \A? + Ned | J 
YA xA, 2 
then - < Apy : ——t— , (69) 
\~ : A? a 
1 
x(w, A w) - 
where A; : — — —— © (70) 
, v“(w,w) (Aw, Aw) 
[hese relations are shown graphically in Fig. 10. 
a, 
a eee 
ARM : } i a 
! i 
' | Xd 
4 ' ! vy 
$— 7 
a 1 6 
S 
asi inky 
Fic. 10. 
Let us define a positive quantity e such that 
| r} \ 0 
Ap a1 when Ap, 
| 
when Apy < 0: 
2(1-+-« 
aiter some reduction. we obtain 
hes 0, A,/« lies in either (—%,a) or (b, +2) (71) 
A | b, (72) 
i 
l le (e7-+-2e) a 
N (73) 
) 1 +-e+,/(e*+ 2e) 
- Lf \ () ) \ ¥ a (74) 
p 
nd \,/x lies in either (—o, —b) or (—a, +0), (75) 
ere a and 6 are defined by (73). 
it is not difticult to show that (72) and (74) are equivalent to the first 


iality given by Southwell (8). 
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15. Case 5 
f(r) (- 


9 
3 ° . fa hl . ‘ 
4 , where « is an arbitrary real constant. The following 
\a A 
method presents another case in which the difficulty encountered j 
section 13 can be removed. In this case the eigenvalues of the transforme 


«= (E45) (i = 1,2,...,n), (7 


matrix f(A) are 


¥ A; 


U 
and they are all positive and equal to or greater than 4. The generalize 


Rayleigh’s principle asserts that if A,, and A, are defined so that 


2 ad min! a i 
K = —-+ —_—_—-+- — 
oi ( X | | xX 7 | = 
A, , «\? (fA; , ~\?) . 
=~ 2 ax! (ft4. 2 
ts) ne | 


A xX 2 A x | 2 
. Pi} . < gq) 70 
the n ( “se Aru 2+5} : (4 
X Ay X W 


((A2-+02)w, (A2-+-021)w) ; 
where Aru = “ 


v?(Aw, Aw) 





These relations are shown graphically in Fig. 11, and we get 











where ; ee. 
a = 148—,/(8?-+28) } 


b 1+-6- (8? 28) fe 
and Aru 4(] | 5)? J 


This is equivalent to the second restrictive inequality given by Soul 


well (8). 
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theorems which provide the bounds of eigenvalues can be derived from 
the generalized Rayleigh’s principle. The next problem would be to find 
the best and the most efficient method for locating the bounds of eigen- 
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16. Final remarks 


It has been shown in the latter part of the present paper that several 


ilues. In the above discussion all the roots on the real axis in the 
\-plane were mapped on the real axis in the x-plane by the transformation 
f(A). The general theory of transforming A-plane into «x-plane would 


e worthy of consideration. 
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HEAT TRANSFER FROM A FLAT PLATE THROUGH 
A TURBULENT BOUNDARY LAYER 


By D. R. DAVIES (University of Sheffield) 
[Received 28 October 1954. Revise received 21 December 1954] 


SUMMARY 

A similarity solution of the heat transfer equation, suggested by Townsend’s recent 
experiments, is obtained for a turbulent boundary layer overlying a flat semi-infinit 
heated plate, provided that the temperature difference between the plate and the fre 
stream is small, and an expression is derived describing the mean rate of heat transfer 
from a section of the plate. The application of this solution requires a knowledg: 
of the observed mean velocity and turbulent shear stress profiles, and is based 
the assumption that heat diffusivity is proportional to momentum diffusivity. T! 
measurements made recently by Townsend of mean velocity and turbulent shear 
stress profiles in flow over a smooth plane wall enable the distribution of temperatur 
in a turbulent boundary layer to be calculated. The theoretical results are compared 
with temperature measurements made by Elias over a heated flat plate, and th: 
comparison indicates a value in the neighbourhood of 1-0 for the ratio of heat 
diffusivity to momentum diffusivity. 


1. Introduction 
RECENT experiments by Townsend (1) on the structure of a turbulent 
boundary layer over a smooth plane wall show that a distance a, may bi 
chosen so that the mean velocity and shear stress profiles may be wel 
expressed as functions of the non-dimensional variable € = y/(2,—2 
where y denotes distance measured perpendicular to the wall and », 
measures distance along the wall from the wind-tunnel entrance. This 
similarity property discovered by Townsend over his range of measure: 
ments suggests that the associated problem of heat transfer from the wall 
to the air stream may be solved by an extension of the method developed 
(2, p. 780) in the case of forced convection in a laminar boundary layer 01 
a plate in a uniform stream. The temperatures of the surface of the flat 
plate, or wall, and free stream are taken to be maintained at constant 
values 7, and 7, respectively, and the temperature difference 7,—, 
taken to be sufficiently small to prevent the appearance of buoyancy effects 
so that the mean velocity distribution is the same as for isothermal flov 
The experimental results giving the similarity distribution over a smoot 
plane wall of mean velocity and turbulent shearing stress enabled Townsett 
(1) to calculate the distribution of the eddy momentum diffusivity or eddy 
viscosity. Following the Reynolds analogy between heat and momenitl! 


(Quart. Journ. Mech. and Applied Math., Vol. VIII, Pt. 3 (1955)] 
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transfer (2, p. 819), we then assume that the dependence of the eddy heat 
diffusivity «,, and the eddy momentum diffusivity « (using Townsend’s 
notation) on the € variable are of the same functional form, and we write 


€x = Ae. (1) 


A formal similarity solution of the ensuing heat diffusion equation for the 
fully turbulent zone of flow is then obtained. An approximate method of 
evaluating the flow in the laminar and transition zones near the surface is 
then used to complete the solution, and formulae describing the distribution 
over the plane boundary of the mean rate of heat transfer are derived. 
A numerical solution based on Townsend’s observed values of € is discussed 
in relation to measurements of temperature obtained by Elias (3) in a 
turbulent boundary layer over a smooth flat plate. Comparison of theoreti- 
cal and experimental temperature distributions indicates a value of A in the 


neighbour hood of 1-0. 


2. Similarity solution of the heat-transfer equation in the fully 
turbulent region of flow over a semi-infinite flat plate 
Let (U’, V) denote components of mean velocity, in a turbulent boundary 
layer, relative to coordinate axes Oz, in the plane of the plate and parallel 
to the direction of the free stream, and Oy at right angles to the plate. 
Then, denoting the mean temperature by 7’, the thermal diffusivity by «, 
ind the eddy heat diffusivity by «,,, the equation describing the distribution 


of temperature as given by Howarth (2, p. 821) is 


eT oT @é[ eT oT , 
l | = —|K —-+€,, —}. (2) 
Cx cy cy coy cy 
The eddy viscosity « in a turbulent boundary layer has been evaluated 


numerically by Townsend from his observations of shearing stress 7, and 
detined by the equationr = «@U/éy. We find that, except for the extremely 
thin region in the neighbourhood of the laminar sub-layer, ¢,, has a con- 
siderably larger numerical value than «, which we therefore ignore in the 


fully turbulent region of flow. We then write our basic diffusion equation 


in the form : - , ann 
oO! ey c eT : 
l SS A —le —I, (3) 
a oy cy\ ey 

which relates to the mean temperature distribution in the fully turbulent 
region 

It was discovered by Townsend (1) that his experimental results for mean 
velocity, turbulent shear stress, and turbulent intensity distributions could 
be expressed as functions of the variable g yx-!, where x measures dis- 


tance from a point 51 em. beyond the entrance to the working section of 
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the wind tunnel. It is convenient, however, to obtain first a similarity | 
solution using the more general form 

€ = y/(ka’), (4 
where the parameters k and q, and a convenient origin for 2, are to be 
determined from experimental results. 

Following the usual Reynolds turbulence technique of taking mean | 
values over a sufficiently long time-interval at a point, the continuity | 
condition for an incompressible fluid applies to the mean flow, i.e. 

oU. av oe 
cael - 0. 5 
ox cy 


We then assume the existence of a stream-function & in the form 


b= x'f(é), 6 
so that U 1 k-lg-Of"(€), i) 
cy 
and V = — 4 ra’ ftqa’téf’, (8) | 
Ca 


where a prime indicates differentiation with respect to €. 
Extending Townsend’s experimental results, we write 


e= KU xPh(é), (9 


where Uj, is the free stream velocity, and A(é) is assumed to be known 


numerically from experimental results. Townsend finds the values 

kK p= I 
to be convenient, but other values are found to be necessary in discussing 
the heat-transfer measurements made by Elias. We assume that similarity 


solutions of equation (3) exist in which 7' is a function of €. We then have 


oT mi eT mr eT 2. osm 
- — — gx jal = -ly-a7” | —— = k-24-727T" ”, 
Cx cy cy" 


€ 


and equation (3) becomes 


iT’ +hT" = — ( bf \ ir" (10) 
KAU 
This equation is easily integrated to give the formal solution 
” kr [ f(é) dé\ | 
[=A | [h(€)]-ex | ; SS 1é + B, (11 
| [m()] °*P)| — EAU, | née) | 
é, é, 


where &, corresponds to the lower limit of the fully turbulent region. 
If the plate and the free stream are kept at constant temperatures T, and 
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larity | 7, respectively, and if the temperature at the lower limit = €, of the fully 


turbulent region is 7'T, the boundary conditions are 
? T->T, asé->0, T->+TYj asé>é,, TT, 


Using the form (11) and these boundary conditions, the temperature distri- 


bution in the fully turbulent region of flow is given formally by 


r | F(é) dé 
— : oe . (13) 
\—TY r 
| P(g) dé 
fs 
& 
em F(s a kr = f(&) d&| 
where (Cc) WC) | exp : - = 
| KAU, . h(é) | 
é 
Si 
lhe parameters k, A, r,q are to be calculated from the experimental results, 
Q nd f(é) obtained from equation (7), when the observed mean velocity 
profile is known 
\n expression for the mean rate of transfer of heat in the fully turbulent 
region per unit width of a plane section perpendicular to the plate may 
w be obtained from equation (13) by evaluating the integral 
Q,(2 pc, | U(T—T) dy (14) 
ie 
) 
ta given station a r,, where p is the air density, c,, is the specific heat 
me | fair, and y, denotes the lower limit of the fully turbulent zone. Using 
LAarity aa | ' , y 7 -_ . 
lownsend’s results, we find that the power law U/U, = 1-72€°} fits the 


| experimental results very closely over almost the whole range of €. Substi- 
tuting this result, together with the temperature formula (13), into the flux 





expression (14 leads to the form 
Q4(x) = 1-72kpe, U,(T*—T)a4 | } | Fdé/[F ag| g015 dé (15) 
p > | S C c 
e é } ¢ 
| oe Ox ox 


noting that the mean rate of heat transfer in the fully turbulent region of 
flow is proportional to a 

This quantity may also be evaluated by an alternative method. The 
upward flux of heat at a point on the lower limit = €, of the fully turbulent 


Zone is cviven by 


oc AecT ley aty Y,, for a given x. 


































330 D. R. DAVIES 





Using equation (13), this expression becomes Howe’ 





K pep AUy(Ti—Th)x?4 | y/ (kat 
ro ’ made 
k | Fdé pons 
: Conse 
bs 
and consequently the rate of heat transfer through a length 2 = 2, of the 
gee 3, Ar 
surface € = €, is given by 
(pre 
at T(q*®__ J! \y1+p-@ 
Q(x) — K pe pal of 1 1 Ty) ; (16 pla 
k(1+p—q) ( F dé Me 
f. (2, p. 


For a consistent diffusion model we see, by comparing expressions (15) and 
(16), that 

l+p—q = 4q, (17) | in the 
and we note that the proportionality constants involved in expressions 
(15) and (16) must also agree. We note from equation (7) that, if U is to 





; ; in th 
be a function of €, then we must write 
q rf. (18 
Nae : ee in th 
Hence, if q is determined from measurements of U/U), the values of p andr 
. . = , the ¢ 
follow from equations (17) and (18): the values of the parameters k and K nT 3 
ee pres , is 
are also to be determined from the U/U, and e/U, measurements. 
‘ . i (20), 
It is also of interest to note here that, if the measurements of ¢ can be , ) 
. cia ci ., | from 
fitted to a simple power law, the temperature distribution may be obtained | 
, aE ae 3 | meas 
in terms of tabulated mathematical functions. Ifh = a&* and U/U, = bé, ) . 
; : : ayel 
then the function F(é) may be written in the form iis 
just 
: k2rb £2+B-x__¢2+B x ; 
F(é) = a-lé-%exp} 4 (5 - —— \. A 
| KadAa(1+ f)(2+B—a)} the | 
and consequently the temperature distribution function may be expressed | ‘oT 
in the form 





T : ° 











whe 
1)—TY fully 
k2rbé2+B8-« -(1+-B) | { k2rb&2+P-% (1+) | 
| KAa(1-+B)(2-+-B—a)’ (2+B—a)| | KAa(1+B)(2+-B—a)’ (2+B—a)! 
_f|__BrbgerPs —(1+8) | 
| KAa(1+-8)(2+B—a)’ (2+B—a)| 
(19 at a 
where the J(x, p) functions are defined and tabulated by Pearson (4). Pr, 


We now require solutions of equation (2) in the laminar and transition 
zones of flow very near the surface in order to complete the problem 
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However, similarity solutions of equation (2) in terms of the variable 
)/(ka®) do not exist, and, moreover, the velocity profile measurements 
made by Townsend and Elias do not extend into the laminar sub-layer. 
Consequently, we adopt an approximate method of procedure. 


3, An approximate method of treatment of the heat-transfer 
problem in the laminar and transition zones of flow over a flat 
plate 
Measurements of the mean velocity distribution very near a plane wall 

2, p. 824) show that the mean velocity forms 


= 


5) at U/U U.yjvu=y, (y, <5) (20) 
7) | in the laminar sub-layer, 
ssions U/U. 3°05+5log.y, (5 <y, < 30) (21) 
: in the transition layer, and 


fe. 5-5-+-2-5 log, yY, (y, = 30) (22) 
oni’ the fully turbulent layer, may be used to give a good representation of 
ey the dependence of the mean velocity of flow on distance from the wall: 
(is the ‘friction velocity’, whose value is determined by fitting the forms 
-_ 20), (21), and (22) to the measured values of U(y) at a given distance x 
sia | from the leading edge of the wall, and v is the kinematic viscosity. The 
; measurements made by Elias and Townsend did not extend into the laminar 
| layer, and we fit equation (22) to the observed values of U at values of y 

just inside the fully turbulent zone. 
| Assuming that the shear stress and the vertical heat flux at a point in 
| the laminar and transition layers are independent of y and equal to 7, and 

‘essed | o Tespectively, we have (2, p. 823) 

IV pe, i Se oll (23) 
where 7 is the ‘friction temperature’. At a point on the lower limit of the 
fully turbulent region this is also given by the present analysis to be 

K pe, AU,( TF -Ty x? 9 
ad Go . « ] (24) 
h ( F dé 
bs 
1" it a given station a v, from the leading edge. 
The temperature at a point in the transition layer is given by the form 
= 2, p. 83] — , 
. : “log, Ao, 1] - Aa} 1 Ho, 
ble 7" A =. 5 











332 





D. R. DAVIES 


where o is the Prandtl number. This becomes 


N-Tt 5 | 
2 : log,(5Ao+-1)-+-5e, (25) 


i A sid 


at y = y,. Eliminating 7. from (23), (24), (25), we obtain 


ee 5KAU,{(1/A)log,(5Ae+1)-+-o% sa | 
1—-T pas oe | , 
kU at-» [ F dé 


é 
Sm 


and finally, using equations (13) and (26), the mean temperature distribu. 
tion is given by 


be T (2 m 
rom _ (qty) (1,73 








Hg) = i 
Th T, ] (7 a 
FP dé) || [ F dé) 4 5KAU,{(1/A)log,(5A0+ 1)+} 
é, é, kU, aq-P ( F dé 
F 
|, BRAT, {(1/A)log,(5Aa-+ 1) +0} 
kULat-» | F dé 


Anexpression for the flux of heat through the region from y = Otoy = 4, 
across unit width of a plane section normal to the flat plate, at a distance 
“% = x,, can be obtained from empirical velocity and temperature forms 
(2, p. 831). We find that the contribution to the flux in the transition zone 
from y, = 5 to y, = 30 is given, at x = a, say, by 


30v/U; 
n 





pC), U(T—T,) dy 


5v/U; 
30v/U> 
a Al JIV\ TY ‘ = = 
pc, (T,— >) U, | | —s05 | 5 log, 


. 
5v/U; 


ly 


Vy 


(5/A)log, {(AoU_ y/5v) + 1—Ao}+-5e ly, (2 
Te ay, (2d 


(5/A)log, (5Ao-+1)+5o0-+kU,a4-? | F dé 


¢ 
S* 


with a similar result for the flux through the laminar region of flow from 
Yy, 0 to 9, = 5. 
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HEAT TRANSFER FROM A FLAT PLATE 


4, A comparison of theoretical and experimental evaluations of 
temperature distribution above and rate of heat transfer from a 
5 heated smooth flat plate into a turbulent boundary layer 

Townsend (1, Fig. 2) shows that U/U, over a smooth plane wall is a 
function of y/x, a being measured from a point 51 em. beyond the entrance 
to the wall, and that e/a) (1, Fig. 10) is also a function of y/x, but he gives 

no temperature measurements with which we can compare theoretical 


results. Elias (3), on the other hand, provides temperature and velocity 


profiles over a smooth flat plate for various free stream velocities and 
si | distances up to 50 em. from the leading edge of the plate, but gives no 
measurements from which e could be calculated. In Townsend’s experi- 
ments the boundary layer is thickened by an obstruction near the entrance, 
nd in the experiments made by Elias the boundary layer is completely 
uminar in the region near the leading edge, but in both cases the flow is 
ver a smooth flat surface and temperature conditions are such that 
buoyancy effects are avoided. This suggests that we may extend the 
ilues of e obtained by Townsend to the conditions of Elias’s experiments 
ind thus enable a quantitative comparison to be made of theoretical and 

measured temperature distributions. 
97 We shall assume that the effect on the flow in the fully turbulent region 
due to the presence of a completely laminar layer near the leading edge of 
flat plate may be allowed for by writing the similarity variable in the 
form y || k(a—a,)“|, where x denotes the distance from the leading edge of the 
plate and x, denotes the distance of the virtual origin of the fully turbulent 

: } region from the leading edge. 

— We first discuss the results associated with the highest free stream velocity 
35 m./sec.) used by Elias, so that any complicating effect due to the 


presence of a completely laminar region near the leading edge is as small 
is possible. The mean velocity measurements made by Elias are rather 
scattered and are not sufficiently numerous to determine exactly the best 
values of x, and g. We find that any value of q between 0-8 and 1-0 could 
be used with a zero x, value to provide a reasonably satisfactory similarity 
variable, and in order to ec ympute the formal expressions (15), (16), (27), (28) 
we choose a value for g of 0-89 with x, = 0. The dependence on x of the 
expression (15) for the flux of heat in the fully turbulent region of flow 
| then agrees with the measured result given by Elias (3); in this connexion 
we also note that a numerical evaluation of expressions (15) and (28) shows 
that 95 per cent. of the total flux at a given x value takes place in the fully 
turbulent region of flow. Consequently the contribution to the down- 
wind flux of heat in the laminar and transition zone does not significantly 


affect the dependence on x of the net result. The numerical value of k is 
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taken to be 1-02, so that the velocity measurements made by Elias for 
¢ = 20, 30, 40, 50 cm. at U, = 35 m./sec. coincide quite well (Fig. 1) with 


the curve, giving the dependence of U/U, on €, plotted by Townsend. 
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Fic. 1. Mean flow distribution U/U,, measured 
by Elias at U, = 35m./sec. 
x 20 em. from the leading edge of the plate 
L 30cm. 
z 40cm. ” 
@ 50 em, 99 ” 
— U/U, = 1-72&" 
€ = y/(1-022°"**), where z is measured from the leading edge of 


the plate. 


The value of the parameter p is now given by equation (17) to be 0-75 
and we consequently take ¢«/(KU, x78) to be a function of the similarity 
variable y/(1-02x°%*). We next determine an appropriate A value by 
comparing with Townsend’s work, in which the similarity property was 

ay), the 
a,). This 


demonstrated by plotting U/U, and e/[(w—a,)U)] against y/(x 
quantities «/U, and y being multiplied by the same factor 1/(x 


suggests that we choose the value of K so that the variable Kx78 coincides 
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ias for | numerically as closely as possible with 1-02x°8® over the range of x con- 
) with } sidered, and we find that the closest numerical agreement for valuesof x in 
id. the range 20 cm. to 50 em. is obtained with K 1-54, the numerical values 
ff Kx®78 and kx®’®® being then within 1 per cent. of each other. Since our 
hoice of k value brings the Elias and Townsend U’/U), curves into coinci- 
lence, we now assume that Townsend’s ¢ curve (1, Fig. 10) can be applied 
in the conditions tested by Elias with ¢/(1-54U,2°7§) and y/(1-02x%89) 
replacing €/|(#—a»)U)| and y/(~—a). The values of h(£) given by Townsend 
ve then used in the calculation of the function F(é) forA = 0-8, 1-0, and 1-2, 


¢) having been obtained by integrating equation (7). 
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Fic. 2. Theoretical and experimental temperature 
distributions. 
v Elias 
| from the leading edge of the plate 
78 — ted @ distribution for A 0-8, 1-0, 1-2 
oe 72 
ul ris measured from the leading edge of 
the plate 
The 6(€) function wasthen ec ymputed from equation (27), usingappropriate 


values of ( and €, calculated by the method described in section 3 and 
This king ¢ 0-72; the calculation was repeated for values of x between 20 cm. 


ind 00 cm., no significant change being obtained in the theoretical 6 values. 
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The results are shown in Fig. 2, together with the measured values of 6 dye | effects 
to Elias (3). The theoretical 6(€) curve, based on A = 1-0, corresponds } transf 
fairly closely to the measured values, and we infer this value of 1-0 for the | the va 
important parameter A, giving the ratio between the heat and momentum 


diffusivities. We then find that the theoretical net rate of transfer of heat | ace 
Q(x), obtained by adding the value calculated for y > y, from equation | The 
(15), with A = 1-0, to the value calculated for y < y,, is in agreement with | 
the experimental result Q/(7,—7)) = 0-56 watts per degree centigrade. | 
obtained by Elias for a length 2 = 50 cm., when Ll), = 35 m./sec.; we als | 
obtain the same calculated value for Q(x) by adding the value computed | 


from the alternative expression (16) for the flux in the fully turbulent regio: i 
of flow, y > y,, to the value computed for y < y,. | 21.1 
Finally, we note that the dependence of the flux of heat on U), is given } 3 F:! 

4 


formally by expression (15) or (16), and in order to evaluate this dependence 
numerically we first need an accurate assessment of the variation of the / 


— 


parameter with ly. When U, = 25 m./sec. or less, the problem is compli- 
cated by the presence of an appreciable completely laminar layer extendin 


from the leading edge of the plate: thus at Uj = 25 m./sec. the velocity | 
measurements for x = 20, 30, 40, 50 em. when plotted against y/(«—2, 


display the similarity property approximately if we choose values of : 
between 5 cm. and 8 cm. and values of q of about 0-89. The points ar 
however, too scattered to determine a reasonably accurate value of x, an’ 
hence of k, and a similar result is obtained for Ul, values less than 25 m./se 
Elias’s measurements are thus not sufficiently accurate to determin 
numerically the dependence of & and hence of Q(a—2,) on U5. 


5. Conclusions 
In this paper it is shown that a suitable similarity solution exists of thi 
heat diffusion equation in the fully turbulent region of flow in the boundar 





layer over a flat surface. Using the mean velocity profile measurements | 
given by Elias for flow over a smooth flat plate and extending to flow ove! 
a plate the mean shearing stress measurements made by Townsend fo! 
flow over a smooth tunnel wall, the theoretical temperature distributio! 
over. and rate of heat transfer from, a section of the plate can be calculated 
We obtain the important result that the best agreement between tli 
theoretical and measured temperature distribution is obtained if we assut 

a value of 1-0 for the ratio, A, of heat diffusivity to momentum diffusivit} 
noting that this result is based on the assumption that Townsend’s obser\ 
tions of shearing stress near a smooth plane wall may be extended to flov 
over a smooth flat plate. This is a reasonable assumption, however, as U 
physical conditions of flow in the two cases are similar, with no buoyan 
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'@due | effects. The numerical value of the expression obtained for rate of heat 
ponds } transfer from a section of the plate is found to be in good agreement with 


or the | the value given by Elias. 


ent 7 
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THE INEXTENSIONAL THEORY FOR THIN FLAT 
PLATES 


By E. H. MANSFIELD} (Royal Aircraft Establishment, Farnborough) 
[Received 7 May 1954] 


SUMMARY 
This paper presents the inextensional theory for thin flat plates. The theory is 
applicable to cantilever and other plates subjected to a normal loading which is 
resisted primarily by the flexural rigidity of the plate rather than by the extensiona 
rigidity. The middle surface of the plate is assumed to be inextensional so that thy 


mode of deformation is a developable surface. 


1. Introduction 

In considering the behaviour of flat plates subject to normal loading it has 
been customary to use the classical small-deflexion theory. It can be shown, 
however, that for thin plates the small-deflexion theory soon breaks dow 
and should be replaced by the von Karman large-deflexion equations (1 
which incorporate the effects of the induced middle surface forces and 


strains. These equations are almost intractable, but for two types of 


problem the equations may be appreciably simplified. 

The first type of problem is characterized by the fact that the plate is 
supported along the boundaries so that the load tends to be taken to a 
increasing extent by the middle surface forces. A typical example of this 
first type is a rectangular plate supported along its edges and subject to 
say, a uniformly distributed load. Such problems may be simplified by 
taking the flexural rigidity equal to zero, though the resulting equations 
first derived by A. Féppl (2), are still difficult to solve. The equations are, 
of course, non-linear because the deflexion is not proportional to the load 

The second type of problem, which is considered in detail here, is charac: 
terized by the fact that the plate acts as a cantilever, being clamped along 
one straight boundary and free along the others. The plate tends t 
resist the applied loading by its flexural rigidity only; and because of tht 
boundary conditions the middle surface forces play an insignificant part 
in resisting the load. A simplifying assumption which can now be made! 
that the middle surface strains are zero, and because of this the simplifit 
theory will be called the inextensional theory since the middle surface 

Acknowledgement is made to the Chief Scientist, Ministry of Supply, for permiss 


to publish this paper. Crown copyright reserved. Reproduced by permission of the ( 
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‘nextensional. One of the differential equations for the deflexion of the plate 


| then reduces to bea) eww 4 


| 


| This equation expresses the fact that the surface is a developable surface; 


a 9 ‘ 9 
Cxcy ox" oy” 


| this follows from the fact that the middle surface is inextensional. The 
ii 

equation does not, however, prove very useful per se as a starting-point for 

solving a particular problem. It is better to start from the fact that the 


surface is a developable surface and to determine the generators of the 
surface from considerations of strain energy. It will be shown that a plate 
oroblem then reduces to the solution of a one-dimensional non-linear 
lifferential equation of the second order. Alternatively, the theory may 
be coupled with an energy technique to give an approximate solution. It 
will be shown too that because the sum of two developable surfaces with 
lifferent generators is not itself a developable surface, the principle of 
superposition for different load systems is not generally valid, but for any 


ne given type of loading the deflexion is proportional to the load. Examples 


| of plate problems to which the ineatensional theory is applicable and which 
annot be classed as cantilevers include: (a) a square plate loaded by 
oncentrated forces at the corners, which deforms cylindrically with the 
generators parallel to a diagonal and does not take up the form of the 
saddle which small-deflexion theory would predict; (b) a strip under end 
moment; it can be readily shown (3) that the anticlastic curvature pre- 
ted by small-deflexion theory has but a fleeting existence in reality; 

his . shallow truncated cone formed from an originally flat sheet bounded 
| radii and ares of concentric circles. This last case differs from the 
preceding ones in that the generators are not parallel but all meet at 
the apex of the cone. When we come to consider the general problem of 
he bending of thin plates it must be realized that if adjacent generators 

| intersect, the sheet bounded by these generators will form part of a conical 


| surface, but the ‘radius’ of the cone and its ‘curvature’ may vary con- 


The locus of the cone centres i the ‘edge 


regress of the developable surface (4). 


uOoUsSLY el the whole plat 


2. List of symbols 


tance along reference axis; a dash denotes differentiation with 
respect to a 
e made by a generator with the x-axis 
nce measured along a generator from an elemental cone apex 
iT the iX1S 


te boundaries 
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M, moment per unit length about a generator 
M, total moment about the «a-generator 
t thickness of plate 


Et?/12(1—v?) 
ce constant introduced in equation (2) 


U strain energy 
F F(a,2,x'), integrand in square brackets in equation (4) 
o outer fibre stress 
F @F/éx, where F is regarded as a function of independent variable 
a, x, x; similarly for other suffixes 
p 72/1 
Bh p—l 


w width of strip 
M applied moment 
P applied load 


3. Outline of the inextensional theory 

The method of solution will now be given to the problem of a thin cant 
lever plate of variable thickness under an arbitrary loading normal to t! 
plane of the plate. The mode of deformation is a developable surface a1 
attention is given to the determination of the directions of the generator 


of this developable surface. An arbitrary distribution for the generator | 
is assumed and the strain energy stored in the plate is then determined | 


The correct distribution is to be determined from the condition that th 
strain energy is a maximum. 

Consider the cantilever plate shown in Fig. 1. Let two adjacent generator 
making angles « and a+da« with the z-axis intersect at the point H so that 
the elemental strip KLL’K’ forms part of a conical surface whose apex i: 
at H. From geometrical considerations the distance along the generat 
from the apex H to the z-axis is given by 


Nr = v'sina. 


3.1. Equilibrium 

An equilibrium condition will now be formed by taking moments abot! 
a generator, for then the component due to the unknown middle surf 
forces vanishes. The moment resisted by the element dy is M, dy andt 
total moment about the «a-generator is therefore 


Ne 


| M,, dy. 


"1 


D, flexural rigidity of plate at distance » along a_ generator, 
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But from the conical properties of the strip the radius of curvature normal 
to the generator is proportional to the distance from the apex, so that 


} 
“ Cc, say, (2) 
where ¢ is to be determined from the equilibrium condition, namely 

! "pb 
M, dy c 7 a | 
‘ J 7 / (3) 


M 


x 


the moment of the applied loads. 
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3.2. Strain energy 
In calculating the strain energy in the plate it is convenient to integrat 


first over elemental strips bounded by generators, and then to integrate | 


with respect to a. Thus 


. 


2 Ne 
U 1 ; M; dnd lf ee " (‘D-In) dn| dx 4 
—1 D n ayn ax - wes § nl) I 
J e 1) ae . 
1 1 ; 
from equations (2) and (3). 

Now 7, and », are known in terms of a, x, and 7,, so that by virtue of 
equation (1) the integrand in square brackets above is a known functior 
of a, x, x’ and the applied loading. For a given loading we can therefor 
write : 

U =} F(a,x,x’) da, say, j 
where the integration is over the whole plate. 

The relation between x and « is to be determined from the condition that 
the strain energy U is a maximum. This condition may be expressed as 
non-linear, second-order differential equation for a: 

Foe fe. — Fe, = ©. (6 


ez’ rr ur r 


The derivation of equation (6) from the condition that expression (5) is 
a maximum is a result of the calculus of variations (5). In equation (6) F 
for example, stands for 0F'/éx, where F is formally regarded as a functior 
of independent variables «, x, x’. Equation (6) is the differential equation 
for determining the distribution of the generators of the developabl 
surface to which the plate deforms, but since the shape and thickness of th 
plate and the applied loading all enter into the equation, no further progres: 
can be made in integrating it unless an actual plate and a definite applied 
loading is considered. Let us consider, therefore, a simple example, ye! 


one apparently that classical small-deflexion theory has not yet solved. 


4. Examples 


4.1. Swept plate of uniform thickness 

Consider a long strip one end of which is built-in at an angle a, as show! 
in Fig. 2. We require the stress distribution when a pure moment is appli 
at the far end. 

If we take moments about a typical generator, we find 


M M sina, 


x 











while f 


Thus f, 
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considerations and equation (1) we have 


x’ sin ¥ 
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we find 


sin A 
since we have assumed D constant) 
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and because this expression does not contain the term 2, it will be found 
that equation (6) may be integrated once to give 





F, = a constant, 
9 9 9 
i‘ w cosec?a , weosec*a) |” . 
whence x 2fy + ——__—_ log! 1 +-- . nn, A... (10) 
| x x | 


The constant K, may be determined from the boundary condition that, 
for large x, « tends to }7 and 2’ tends to infinity, so that on taking the limit 
K, = w'*. (11 
The further integration of equation (10) is best done graphically. If we 
introduce a parameter p equal to »,/7, we can obtain a parametric relation. 
I peq 12/1 I 
ship between a and 2’, namely 
a’ I 


w  vplogp | 


vege 
p—l 


and x will then be given by the relation 


and x sin-! 





x ( x’ da. 


Once the (a,2) relationship has been determined the generators of thi 
deflected surface may be drawn, as, for instance, in Fig. 3, where they hav 
been drawn at 5° intervals. 


4.1.1. Stress distribution in the plate. In general, once the distribution oi 
generators over the plate has been found the stresses may be found readil\ 
from equations (1), (2), (3) and from the fact that 
6M, 

t2 


oO (13 

The stresses so found are principal stresses and will be directed at right 
angles to the generators. The stresses in the direction of the generator 
will be v times those given by equation (13), owing to the fact that th 
generator lines remain straight. 

For the particular problem of the swept plate under a pure moment thi 
outer fibre stresses some way from the root assume a constant direction al 
a magnitude given by 6M 
wt? ’ 
and it is convenient to regard this stress as our unit, so that the stress 
elsewhere in the plate will be numerically the same as a stress concentratl 
factor, and can therefore be more readily appreciated. 

















so t! 


and 
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fou 70° 1s 80° 85° 





\ | 
| 


Width of plate w 








From equation 3) we have. on simplifying, 
M sin « 
Diog p’ 


so that from « quation 


M w sin « 








M uw ” log Pp 
= wsin a ; . (15) 
cv sina+(n n,)}log p 


ressions occurring here are known from geometrical or other 
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considerations. The maximum stress in the plate occurs at the junction of 
the trailing edge and the built-in edge (where « = a, and 7 = »,). Taking 
o,, as unity, it will be found that the maximum stress is given by 


o = Vp, 


so that from equation (12) the maximum stress is related to the building-in 
angle a, by the formula: 


=< gin-! 
% = sin 





20 log "| } 


a (16 
o*—1 


Some values of a, are given in the table below: 











[t will be seen that as a) decreases o increases rapidly, but as will be seen 
from Fig. 3, where contours of o over the whole plate have been drawn, 
these high stress values tend to be very localized; for example, when ay is 
55° the peak stress value of 5 has dropped to 3 in a distance of 0-035w 
This type of phenomenon is, of course, characteristic of stress concentration 
problems in general. 

In considering Fig. 3 it is important to note that since a generator remains 
straight it may be regarded as clamped so that any one of them may be 
regarded as a built-in root. Thus we have obtained the solution for the 
swept plate built-in at an angle a). For example, if a, is actually 45°, our 
solution is represented by the plate to the right of the 45°-generator in Fig. 3; 
we simply ignore the additional information contained in the rest of the 
figure. The fact that a whole class of problems of cantilevered plates can be 
solved in a single analysis is an important and useful characteristic of the 
inextensional theory. 


4.2. Swept plate of varying thickness 

The plate considered in the example above was of constant thickness, but 
the analysis is not necessarily over-complicated when the plate thickness 
varies. We can illustrate this by outlining the analysis for the same prob- 
lem except that a cross-section of the plate is now diamond-shaped, being 
of thickness ¢ at the mid-chord and tapering uniformly to zero at the leading 
and trailing edges. 

Equations (7) and (8) will be the same since they depend on the loading 
and geometry. Along a generator the thickness varies linearly from zero a 
7 = 7, tot at 7 = 4(y, +7.) and then linearly to zero at 7 = yj». Now the 
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fexural rigidity varies as f°, so that along a generator the flexural rigidity 


arieS as follows: 


3 
' a 
1O1 4 ” >(14 > No) we have spi ] | 
2-1 n 
In (17) 
3 
Ne 7 
101 (1) Ne) ] No we have sp! 12 1 | 
a 
YD | | | 
term dy occurring in equations (4) and (9) therefore becomes 
72 
8D j 3 5 — 
11) ie ad (> 1)’ di 
3 / 7 
1) F 1) J n 
Ami +7 


reduces to 





SD : 2No 4 +o 5D(n, +n. 
u log| ; ni log| neta) _ Wh 7) (18) 
mi)?” 1 2 211 No—-N1 
nd the function F is finally determined by the equation 
VU 20 +2’ 3 C+ a’ ee 
(2¢-+2 )* log : x "og ; 5O*(C-+a) (19) 
SDF sin4a (+2 x 
\ Ss Ly 
r¢ ( Sw cosec*x. 


lation (1 


does not contain the term 2, so that, as in the previous case, 


the differential equation for the generators integrates once to become 


F a constant, 
20 o( 2 vr’) +2" log 2§2724 Ifa’ +-2’*\log(C xv’) —3l? 
K, sinta x”, (20) 
nd, for the same boundary condition as before, we find, after some algebra, 
3% 
K, , (21) 
. 3w4 
ugh equation (20) is fairly involved, the relation between 2’ and ¢ (and 


in obviously be found by graphical and numerical methods, 
e fir calculation for the stresses then follows a similar course to 


of the first example 


ng | 4.3. Triangular cantilever plate under a tip load 





lhe plate with a typical generator is shown in Fig. 2. This example has 


een chosen because the resulting differential equation for the generators 


the se | order in x and «, and a different technique is needed for 
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Taking moments about a typical generator gives 


M, = Pxsina, 99 


(<2 


while from geometrical considerations 


Hy = 2’ sina 
. 9 
xsind /> 23 
N= M4 a 


“ ln sin(a+¢) J 





so that for a plate of uniform thickness we have 


F P2z? sin?a 
: — ——_—_—_——_— —— (24 
xsind 
Dlog! } + —-_______ | 
| © 2’ sinasin(a+¢)| 
Substituting equation (24) in equation (6), we have an equation of th 
second order in x and a, but it may be written in the form 


dp cos ¢ + cos 28 


N sin 28 % 
i 
‘| sind e 


l : 
Oe siy 


where - 4u2(1-4 p)log(1 | i) 
7 (2+ p)log(1+p.)—2p 

(1+ p)log(1-+ u)l2 + al MoI | n)| 
( | Mm 
2 = ——— -—_——— —_—__— ————$______—_— 
(2+-pn)log(1+p)—2p 


2x sind Ne— 1 
x'(cos d + cos 28) "1 


al 


B = a—}(7—4). 28 


The advantage of introducing the function » is that the differentia 


equation relating » and £ is of the first order. The angle f has been intro- 


duced so that the boundary condition at the tip reduces to 
B=0 at gp = 0. 29 


The method of integration of equation (25) subject to condition (29) wil 
be given shortly. It will be seen that this primary integration is the onl 
real difficulty since, when the relationship between , and f is determined 
the formal solution for 2 and 8 is straightforward: from equation (27 
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x 2sind r 


z (cos d + cos 2B) j 
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which may be integrated to give 


(f  _2singdp | 


r Cexp mm Yi (31) 
I. (cos b + | cos 28)| 
0 
» method of integration of equation (25) is a combination of exact 


I ds und numerical step-by-step and trial-and-error methods. 
When n the limit is taken of equation (25) as » and 8 tend to zero, it is 


und 1Q ” TS 
(“") sind (74 6 cot? 34)| 
dy ; | 24 | 


Say. (32) 


\ step-by-step integration of equation (25) now proceeds as follows. 


\ small value for « is chosen, say 0-1, and a corresponding approximate 


ie of i.e. O-ly) determined from — (32). Substitution in 
uation (25) yields (df 1. The value of 8 at » = 0-1 may then be 
ised by trial and error until 
0-1 dp 
me Oy (5) | (33) 
2 | dy oil 
t | 0-2 may be determined approximately from the 
| 0-1(/dp dp 
~ (Boa ipa + (2 \ (34) 
2 \\delo1 dy. o-2) 
hich (dB/d1)y.. is obtains ee d)o 
I8du)),. A revised value of (d8/du).. may now be obtained from 
ition (25) and thence a revised value of (f),.. from equation (34). Further 
ues are found in a similar manner. 
When the above calculations are completed there is a known relationship 
tween nd ind from equation (30) we can regard 2’/x as a known 


Now the expressions in equation (30) tend to infinity as B 


ds to ind it is necessary to introduce 7'(8) such that 
ae 24 tan id 2sind — 
T'(B - fs (35 
3 (cos d +- cos 28) 
y » I is 6 tends to zero. 
E juatiol 50) now hece mes 


2ytanid ,, — 
= T'(B), (36) 
B 


which may be integrated in the usual manner. 


€ limiting value 7'(0) cannot be obtained numerically from equation 
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(35) and must be found as follows. Differentiating equation (25) and taking 
the limit as 8 tends to zero, it will be found that 


]2 
(' *) =, (37) 
dp? 0 


and taking the limit of equation (35) gives 
‘mm l z 

T(0) = —- ( tan}¢é = tan 4¢. (38 
dp? 0 : ‘i 


Equation (37) may be used to give a more refined version of equation 


(33), i.e. (8) _ Q-] 2 dp 0-] ata | 3 
"= - \" Al , (idea 2 dy.” o) - 


The distribution of generators and principal stresses is shown in Fig. 4 
for the case where ¢ is 60°. The solution obtained applies to the complete 
range of sweep-back angles w. 


4.4. Triangular cantilever plate under a uniformly distributed load 

The analysis for this case is almost identical with the previous case except 
that in equation (25) the expression —2Q becomes —6Q. Solutions for 
triangular plates are given in greater detail in (6) and (7). 


5. Verification of the theory 

It would be ideal if the results of this inextensional theory could be con- 
pared with a number of exact solutions of cantilever plate problems. The 
only such problem for which an exact solution is known is the rather trivial 
case of a long strip under pure moment. For this problem it is found that 
the inextensional theory is appreciably more accurate than small-deflexion 


theory if the loading is such that 


("") o 3 (40 
t} E 


and as the loading is further increased the results of the ineatensional theory 
approach asymptotically those of exact theory. To fix ideas as to th 
numerical magnitude of the parameter in equation (40) it will be noted that 
if w/t is 100 and o/F is 0-004, 


which is well in excess of 3. 
As a practical check on the adequacy of the inextensional theory a series 
of experiments was carried out and in all cases agreement with the theory 


was good. The first experiment was a check on the tip deflexions of six 








right- 


hypot 
lthou 
dicted 
the dis 
clamp 
compr 
remov 
strain. 
and sl] 
except 
and 8! 
ment | 
plates 











THE INEXTENSIONAI 
king 
ontours of con 
s | 
| 
ple 
‘1 
x 
x 
T 
8 Gene 
LO" 
g ingied triangulat plate 
t Sé In all cases 
ug) e deflexions wert 
ustril of the gene 
ped g its hypote1 IS¢ 
eSsl rrace painted 
ed ippeared 
\ ing to the ine 
sho cide with é 
ion of the 2 
sé S ion conditi 
ads ipplied at 
es Wi fferent t/w ratios 











tors for a triangular plate of angles 90 








« THEORY FOR THIN FLAT PLATES 


he 





351 


(e) 


0 


M97 
at: 1 
stant = tan4 p 











\ 
aN 
\ ‘ 
0-9 \ \ \\ 
\ ie, \ \ | 
\ a 
\ \ rm | 
\ \ \ 
_ —_. ™ 
—_ 4\ ~ 
) X \ -=S 
» A 
/ \~/ \ ae” I-7 
rat \ / \ X\y - 2:0 
orors \ \ /\ 
Y \/VYe25 
- Mb Ade 3.0 
Fic. 4. 


‘s, clamped at different angles along the 
load-deflexion relation was sensibly linear, 
ibout 1-2 per cent. greater than those pre- 
ry. The second experiment was a check on 
40 .20 
The plate was loaded at the tip and the 


with a strain lacquer. When the load was 


ie strain lacquer along lines of principal 
nsional theory these lines should be straight 
generators. Agreement was generally good 
’ corner, where the stresses would be small 
s would tend to prevail. In the third experi- 

tips of five cantilever equilateral triangular 
The smallest value of t/w tested was 0-0025 
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and agreement with the inextensional theory was apparently exact over the 
complete range of loading. With higher values of t/w it was possible to 
detect an initial decrease of about 20 per cent. in the stiffness of the plate 
under small loads. As the loading increased the stiffness rapidly assumed 
the value predicted by the ineatensional theory. In all these experiments 
the plates were cantilevers and the difference between the ineztensional 
theory and small-deflexion theory would be small. 


6. Conclusions 

A simplified large-deflexion theory has been presented for thin plates 
subjected to normal loading. This inextensional theory is applicable t 
plates in which the loading is resisted primarily by the flexural rigidity of 
the plate. In the theory the solution of a plate problem reduces to the 
solution of a one-dimensional non-linear differential equation of the second 
order. Alternatively, the theory may be coupled with an energy technique 
to give an approximate solution. The principle of superposition is not valid 
for different load systems, although for any given type of loading th 
deflexion is proportional to the load. For thin cantilever plates under very 
low loads the difference between results obtained from small-deflexio: 
theory and from the inextensional theory is not large; for higher loads the 
results obtained from the inextensional theory approach asymptoticall 
those obtained from large-deflexion theory. 
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THE EQUIVALENCE OF CONTINUOUS AND DISCRETE 
MASS DISTRIBUTIONS IN CERTAIN VIBRATION 
PROBLEMS 


By KR. K. LIVESLEY 


(Computing Machine Laboratory, University of Manchester) 
ived 21 October 1954 


SUMMARY 





Several numerical methods for determining the natural frequencies of non-uniform 


| ifts in torsional and transverse vibration require the replacement of any con- 
nuous distribution of matter by an ‘equivalent system’, in which the mass or inertia 
neentrated at particular points. This paper discusses the error introduced by 

s process by considering a few simple cases, in which both the continuous and 
liserete systems can be treated exactly. Following a brief statement of the 

s obta | by previous authors, an expression is derived for the frequency 

rin th e of a uniform simply supported beam in transverse vibration. It is 

that if the continuous mass distribution of the beam is replaced by a set of 
entrated masses, the errors in the natural frequencies are propor 

erse fourth 7 ver of the number of segments into which the beam 


1. Introduction 
ENGINEERING design work often requires an investigation of the transverse 
torsional vibrations of a straight non-uniform beam or shaft. In most 
ses attention is concentrated on the values of the natural frequencies, 
he precise nature of the associated modes being of less importance. In a 


im turbine, for instance, it is essential that the critical frequencies of 





the rotor asse mbly should be outside the normal range of running speeds, 
t the actu 


| modes of the vibrations which occur are rarely of interest. 
The various schemes of numerical analysis which are available may be 
issified according to the way in which the displacement of the shaft is 
ed. In the first group the displacement is represented by a linear 
ndination of certain independent arbitrary functions, the latter being 
sen so as to satisfy the boundary conditions. The methods of Ritz 


1, p. 370) and Galerkin (2) are examples of this approach. In the second 


roup ti 


: e problem is reduced to one with a finite number of degrees of 

ireedom by the introduction of an ‘equivalent system’ of concentrated 
isses and inertias. The motion of the shaft is then treated entirely in 

terms of the displacements of these mass points. This idealization is required 

(Quart. Journ. Mech. and Applied Math., Vol. VIII, Pt. 3 (1955)] 
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in the method of analysing torsional vibrations named after Holzer (j 
p. 261), and also in the extension devised by Prohl (3) for treating th, 
transverse case. 

The methods of Holzer and Prohl are both ‘exact’ in the sense that thy 
give all the modes and frequencies of the equivalent system to any requin 
degree of accuracy. These frequencies, however, will not in general | 


identical with those of the original continuous system, and differences ma 


also occur in the actual modes of vibration, although these are usually Jes 


important. In most cases only the lower frequencies are of interest, so th 
the problem facing the practical computer is the determination of t) 
simplest equivalent system which will give a certain number of the low 
frequencies to some specified degree of accuracy. 


It is apparent that a quantitative expression for the frequency erro 


mentioned above can only be obtained in cases where both the continuow 


system and its discrete equivalent can be analysed exactly, and su 


solutions are difficult to obtain unless the shaft is in fact a uniform on 


carrying a uniform mass distribution. Results for simple systems of this 


type are of little value in themselves, but they may be used as a roug 
guide in estimating the errors in more complicated cases. A few of thes 


solutions are discussed in this paper. 


2. Systems obeying the wave equation 

Perhaps the simplest case in which both solutions can easily be found 
the uniform elastic string in transverse vibration. This problem was ¢ 
pletely solved by Lagrange, and an account of the analysis may be fow 
in Rayleigh’s Theory of Sound (4, vol. 1, p. 172). He considers an elast 
string of line density p, held by a tensile force 7’ between two fixed point 
distant L apart. He reduces the system to one with a finite number of ¢ 
grees of freedom by dividing the string into m sections and splitting the mas 
of each section into two equal concentrated masses at the ends as show 
in Fig. 1, thereby obtaining a set of m—1 moving mass-points each 


mass pl im. 





c } 4 
e—L/m—> 
pLl/em 





ria. I 


If ,,w, is the natural frequency of.the sth mode of vibration, Rayleigh show 


that , - 
x (2m L)sin(s7 2m) | (7 p), 


and hence that 


5) ,/ 1 W, sin(sa/ 2m) (sz 2m), 
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here ,.w, represents the true natural frequency of the continuous system. 
ng t Thus 11 s s small 
Ww. Pe (77/24)(s/m)*. (3) 
In this particular case an s shows that although the frequencies are 1n 
ror. the rmal modes f the discrete case are in fact correct. 1.e. if Y, 


represents the modal displacements of the mass points. then for the sth 





. ) {.sin(ars mm) (i = a m—1) 
Similat ysis, applied to the torsional vibrations of a uniform heavy 
ft. has recently been given by Duncan (5), who also considers the more 
seneral problem of the non-uniform shaft. In the latter case he shows that 
rror in ¢ frequency remains inversely proportional to the square of 


4 


number of sections, but he does not discuss its magnitude. 


| [tis worth noting that D in uses a different method of concentrating 

nerti le replaces e distributed inertia of each segment by an 

| amc placed at tl tre. Fig. 2 shows how this scheme compares 

the met lused by R eigh. In the case of a uniform string or shaft 

l be f i that the 1 iethods give identical results. It is clear that 

oh yroach has some advantages when dealing with shafts formed 

n sections of different diameter, as in such cases concentra 

it the ends of the sections greatly simplifies the calcu 
10 nal ric es 
ieee eee 


+-—Z,/m ——»| p/m 


occur in any finite-difference method of 


ues of a partial differential equation. The wave equation 








endent sp riables, associated with the problem of the 

ting rane, ha treated from this standpoint by Milne (6 

237). He scusses recta lar regions only, obtaining results analogous 
one é ibove 


3, Systems obeying the vibrating beam equation 


V 5 considers the problem of a uniform cantilever in trans 
vibrat nd shows from numerical calculations that the inverse 
requency error still applies. Although it does not appeat 

e€ to of Lan exact ilvtical solution to the discrete problem in 

se, § solution can be found when the beam is simply supported. 
ethod used in the following analysis is very similar to that emploved 


Rayleig] is treatment of the vibrating string. 
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[f we consider a beam of length L, line density p, and flexural rigidity 5] 
simply supported at its ends, the natural frequencies of vibration w, ay 


iven by 9.9 9 7 
given by w, = (7s?) L?),/(HIg/p) (s Sy ae 3 4 


and the corresponding normal modes are 

Y (x) = A,sin(szx L). 
If we now divide the beam intc m equal sections, concentrating the mags oj 
each section in accordance with Fig. 1, our mass distribution becomes oy, 
of m—1 equal masses p, spaced at equal intervals a, where xp = pl/m an 
a L/m. We let y, denote the displacements of these mass points, wher 
r runs from 1 to m—1, as shown in Fig. 3. 


“2. 

| “Se 
= 
m-2 m-| 7m m+\ m+2 








Our aim is to find a solution y,(t) which satisfies exactly the equation 
of motion of the mass points and the boundary conditions at 7 = 0,7 
These conditions are most easily dealt with if we treat the beam as part 
an infinite beam of the same material, as shown in Fig. 3. By introduciny 
‘image points’ in this way, it becomes apparent that the conditions of z 
deflexion and bending moment at r = 0, r = m, are satisfied exactly | 
any set of displacements y, for which 

Y¥_, —Y;, Yoam-r> 
where 7 may now take any value. 

In order to write down the equations of motion of the masses we requ! 
the following result, which the present author has not been able to fi 
explicitly stated elsewhere, although a similar idea is mentioned | 
Southwell (7, p. 183). 

LemMA. [fa uniform beam of flexural rigidity El is loaded by concentrat 


loads W),.....W,, the load spacing a being constant, then for 2 <r <n-? 


r 


(Y,-2—4Y,- + 6y,—4Y, 3 +Yr42) = W_1+4W4W,.. 
a 

Consider the section of beam shown in Fig. 4, where F and M are' 
shear and bending moment acting on the beam just to the right of' 
load W,_,. For 0 < a < 4a, the deflexion of the beam satisfies the equat 
dy 


da2 


El M + Far+-W,_,{x—a}+ Wi{a—2a}+ W,{2—3a}, 


r 
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‘\ being used to indicate terms which are only considered when 


if0O<z<a 


4 a 0) 
, i 


’ 


if a Yr. 














a 
Yr+t Gree 
Wr 7 Wre2 
Fic. 4 
O S) we »b 
Va Fa if WW ’ 
i ay "tx—2a'3 
6 6 6 
W.., ,, Aa? Ba? 
yal sa} x (9) 
6 6 6 
B are arbit constants of integration. Substituting the 
ilues of x in (9) we obtain 
18M /a+-64F+27W,_,+8W,+W,,,+44+B, 
27Mja+27F+ 8W_,+ W 34+ B. 
y 12M a SF Ws 2A +- B, (10) 
“ 3M /a F A+B, 
B. 


fourth difference of both sides of equations (10) we obtain 
result (7). It should be noted that this equation represents 
on of the original differential equation. Replacing the right 


6H g1Ves 


) WwW 


(YY, > 47/ , « 
- . (11) 


6Y, 4 i101 Se 
a’ a 


El 
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which is of course the normal finite-difference form of the differential 
equation for the continuous-load case, with the error term omitted. 


3.1. Using equation (7) we can write down immediately the equations 
of motion governing the mass points. It is important to realize that we are 
comparing exact solutions of two different physical systems, so that the 
approximate expression (11) must not be used. The equations of motior 
are 

6EL ; ‘ 5 

aa (%,-2—4Y,_1 + 6Y,— 44-41 TY p49) = q ue tee t+ Y ria), (12 
where may take any value, provided that the y, satisfy equation (6 
Putting y, )_e’@t. we obtain 
peaPa* 


) 2 1)" 1 6Y, +Y,. 1 2 6Elg ri r : 1 


Following Rayleigh, we try a function ¥, giving exact agreement witl 


the continuous case at the mass-points, i.e. ¥, = Asinrf. This clear 
satisfies equation (6) provided that mf = sz. where s is a positive integer 
para” 


Substituting in (13) and putting & , we obtain 


6Elg 


4-k)sin(7—1)68 


sin(r—2)6 


(6—4k)sin rB—(4--k)sin(7+- 1)8-+sin(r+2)p Q), 
and hence 


2 


2 sin vB cos 28 — 2(4+-k)sin r8 cos B— (6— 4h)sin vB 0, 
since by hypothesis sin rf (), 
? cos 28B— 2(4-+-k)cos B-+-6—4k 0 
16sin4(p 2) 
whence h - — 
6—4sin*(f 2) 
Substituting for /. 8. uw, a, we finally obtain 
m= 4sin?(sa7 2m . ~ 
Ww Ray) 7 ViALg Pp); \+2 


L? \ | 2 sin?(s7 2m 
and comparing this with equation (4) gives 
| 


1 —2 sin?(sz/2m)}- 


sin(sz, 2m) 
(S87 
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pecomes 

~ J (sa 2m) 0-0676(s m)?, (16) 
Duncan's method of mass concentration 

result. We see that not only is the error an 


ut that it is considerably smaller than in the 
i, 


2, 8 we 


For instance, if we take m 


{ 


for the lowest natural frequency of 


n 


6\ /(Elg 9-859 /(Elg 

(5 ),/| p Sis } 

ict solution of the continuous case, 
Elgq/p) 9870 (Ely p), 


0-1 per cent., as against 4-5 per cent. in the 


o Case. 
results for the uniform cantilever we can 

mship 
~ 0-2(s m)*, (17) 


rote 


“| by Dunean) 


accuracy corresponds to an error function 


of 15 segments per com- 


(he surprisingly small error in the case of 
learly due to the fact that the modes for the 
those of the continuous system, while the 


satisfied exactly. It is unlikely that such 


er types of support condition, so that in 


expect an error function similar to that 


1 cases where the modes agree at the mass 


cree elsewhere, and the stress distributions 


the Two Cases. 


In the numerical example 


MN parison between the two modes shows 


I 


lacement at the centre is 0-4 per cent. 


7} 
a 


moment at the same point is 5 per cent. too 


{ 


mass which does not have a great effect 


produce a relatively large change in the 


IO PY 
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THE PRINCIPAL FREQUENCIES OF VIBRATING 
SYSTEMS WITH ELLIPTIC BOUNDARIES 


By S. D. DAYMOND 


(Department of Applied Mathematics, University of Liverpool) 


teceived 17 June 1954] 


SUMMARY 


In the equation V2d(a’, y D(x, y) 0, which occurs in the theory of vibrations, 
quen f the normal modes are, with a convenient choice of units, equal to 
square roots of the eigenvalues vy; appropriate to the boundary condition. Two 
blems are dealt with here in connexion with elliptical boundaries, all of which 
se a constant area 7. In the first problem the boundary condition is db O and 
ely quency Is 6 I he second the boundary condition is Ob on 0 and 
princi} frequency is yu. The numerical values of .A and vy are determined 
0(0-1)1 and for (1 ' 1(0-1)0. By applying a simple theorem of exten- 

is then possible to find such values for an ellipse of any area and eccentricity. 

lhe method is an inverse one in which zeros (€)) of the modified Mathieu func- 


Ce,(€,q) and of (d/d&)Ce,(€,q) are evaluated, when the parameter g takes each 


fa conveniently chosen increasing sequence of positive values. A real zero of the 

runct exists for any ich q; but there are no real zeros of the latter if 
ORG 

It is found that as e increases from zero to unity vA increases from 2-4048... to 


ty, insuch a way that A, (1—e*?) > 77/4, while sp decreases from 1-8411... to zero. 


Sections | and 2 contain, respectively, brief accounts of the relevant properties of 


f 


genvalues and of the theory inderlying the solutions of the two problems. Values 


found in section 3; in addition, brief references are made in this section to 
ipplicat fan asymptotic formula, and to the accuracy of a known approxi- 
inal formula for VA. Values of vu are determined in section 4, Finally, 


tion 51 brief note on the €-derivative of the function Ce,. 


‘he notation used for the Mathieu functions is in general that to be found in (1). 


sional ve use that accompanying the tables (2); thus, we frequently write s 
tq and we adopt in secti 5 the method of normalization explained in (2) as 
8 simplhif he use of these tables in evaluating the derivative mentioned above. 


|. Some elementary properties of eigenvalues 
LET d(x, y) sat isfy the equat ion 
V°¢- vd 0 (1.1) 
ina region R bounded by a closed curve C, where x,y are cartesian co- 
rdinates and v is a parameter. Two problems, arising from two different 
boundary conditions, will be considered; they are: 
mn : 
Prosiem I, in which d 0 on C; the eigenvalues {A,' are such that 
U Ao A, ‘ A, 


{Quart. Journ. Mech. and Applied Math., Vol. VIII, Pt. 3 (1955)] 
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and ProBLeM [I, in which ed/éen 0 on C; the eigenvalues {y;} are suc 
that 
O = py < Hy < py 


The property A; > yp; is well known. Moreover, proofs of the interesting 
property A, > , are given in (3) and in the appendix, by G. Pélya, to that 
paper. In the theory of wave-guides (1, ch. xviii) the significance of 
the latter theorem is that for an arbitrary but constant cross-sectioy 
the absolute (i.e. longest) cut-off wavelength is 27/\u,, and hence that the 
dominant mode is an H-mode and not an E-mode. 

Again we have the further properties (3) that 


(i) if Rand R# are the areas of similar regions, corresponding eigenvalues 


of any order are such that 
A, /A; = fp; R/R: (1.2 


(ii) of all regions R enclosing a prescribed constant area the boundin, 


curve C for which A, is smallest, and also yp, is greatest, is a circk 


The truth of (i) follows at once from (1.1) by changing the variables to Z,i 
which are respectively proportional to x, y. 

In (ii) let the magnitude of the constant area be z. For the circle, A, 
2-4048... x, Say, and vy, is 18411... 8, say, where these numbers ar 
respectively the smallest positive zeros of J,(7) and J}(x), in the usua 


notation for Bessel functions. 


For a rectangle whose sides are in the ratio k:1 (k 1). we have 
ro am(ktk-1) > 27 
and by ah 77. 
Thus, as k > 0, A, © in such a way that A,k — z, and p, > 0 like &. 


Consider now a family EF of ellipses of which each member encloses a! 
), Con 


9 


area 7. When the eccentricity is ¢, the ratio of the axes is ,/(1—¢ 
sideration of the above cases of the rectangle leads usto expect that A, , (1 


9 


varies from a? to a finite limit J, and that y.,/,(1—e?) varies from /* 


+ 


a finite limit m™,, 
that /, 7/4 and m, 0. 
It is clear that the eigenvalues corresponding to any given ellipti 


boundary # can be found from those appropriate to the similar memb 


of E by applying (1.2). In the subsequent treatment of the problems 0! 


evaluating vA, and vz, for such a family of ellipses, over the ranges 0! 
and ,(1—e?) stated in the summary, we shall dispense with the suffixes 
so that VA and vy will represent the principal frequencies to be determin 
in the Problems I and I. 


as e increases from zero to unity. In fact it is shown later 
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2. Basic theory and method of solution of the two problems 


Let the region #& be bounded by an ellipse of eccentricity e and area 7. 


The semi-axes are (1—e? We express (1.1) in elliptic coordinates given 
by Ly c cosh(€+-in), 
where Cc? |—e? (l1—e* e*(1—e*)-*. Writing, in the usual way, 
X(é)} it follows that X and Y must satisfy the Mathieu equations 
\ a—2q cosh 2€)X 0), (2.1) 
) a 24 cos 2y)¥ 0), (2.2) 
where the “characteristic number’ a is a separation constant and 
tq . ve" ve*(1—e?)-:. (2.3) 
the boundary we have 
cosh €, ( c | /cosh €, sinh &). (2.4) 
Now let i jexp &, v vq exp(— &€o). 
nd hence “ue vs cosh &p. (2.5) 
It follows s/e v( 1—e?) (2.6) 
ha —s v( 1 —e?) (2.7) 
| ") f (2.2) are the Mathieu functions of the first 
f det 
se (7). Y) (an SO. 8, Z....k. 
those 2.1) are the odified functions 
C = ce, (2 and Se,,.,(€.4q) USC 4 (0&.q), 
vhich are pectivelv evel nd odd functions? of €. 
he bo ry condition of Problem I requires that 
C QO and Se,,.,(&. 7) 0 


we require further that the appropriate zero in q 


lead ilue of vA which tends to a 2-4048... as &; > 0. Now 
modified function only, tends to J,(vA), the Bessel 
t101 r zero, as ©. Thus. on making the substitutions 
eX] exp nd gq A/2 sinh 2€, Ap?/(p*—1), 
A pre ~ 

UM 7 cosh 2 /-cosh 4€ ie9(4 cosh 68 7 cosh 2€) 
(2.8) 
2 xdd modified functions are respectively Je, s, €) 


roduction. 
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can be reduced to the form 


x 


Jo(vA)+ 2, (p* 1)-*f,.( VA), 


where the /;, are power series. Although these series are only slowly con 
vergent, even when vA assumes its smallest value a, and therefore of 1 
use for computing purposes, the above form makes it clear that a direct 
solution of this problem would require the evaluation of parametric zerms 
of the function (2.8), for prescribed values of &. 

The boundary condition of Problem II requires that the &-derivatiyes 
of the functions Ce,, and Se,,,, should be zero when & = &). Also the 
corresponding parametric zero should lead to a value of \j which tends 
to 8 = 1-841]... from below (see property (ii) of the eigenvalues) as &, + x 
Now (1, p. 368) both (d/d&)Ce, and (d/dé)Se, tend to J{(vj) as € = € 


Closer examination on the lines of the preceding paragraph, however 


> 7 


shows that the frequencies which depend on parametric zeros of (d/dé)Se, 
begin to increase as the eccentricity increases from zero, while those whic! 
depend on zeros of (d/d&)Ce, possess the required property that vj for th 
circle is the largest. Thus a direct solution of Problem IL would involv 
parametric zeros of the function 

| (d/dé)Ce, |e_¢, = sinh €,—3q sinh 3€, + J,g?(—3 sinh 3€, +3 sinh 5é,) 

E=§ - 

The method of solution 

In either problem the direct determination of parametric zeros for giver 
values of €, is impracticable, and the inverse procedure is adopted in whic! 
zeros are determined for prescribed values of s (= 4q), as below. 

The modified Mathieu functions Ce, and Ce, are expressible (1, pp. 159 
160) in series of Bessel functions of argument w — vscosh€&. The coefii- 
cients,t functions of s only, are extensively tabulated in (2). Thus, with 
the aid of (2) and (4), we evaluate for each one of a sequence of values of s 
the zeros wy = vscosh €, of Ce, and of (d/dé)Ce,. The values of e and of: 
(i.e. A in Problem I, u in Problem IT) which correspond then follow from 
(2.6) or (2.7). These corresponding pairs of values of e and vA, or vp, are 
of course at unequal intervals of ¢. Four-decimal values for the regula! 
intervals described in the summary are derived from the above values by 
the method of linear cross means.{ An account of this method is to be 
found in (5), pp. 84-86. 

+ The absolute values of the coefficients depend on the method of normalization of th 


ue of (1) and the Dy, . of (2), where p 0 or 1 here, are propor 
tional for a given value of s. The factor of proportionality is unimportant in the calculatior 


functions ; however, the A 


of zeros. 
¢t The author is indebted to Dr. C. W. Jones for pointing out the suitability of ths 
method. 
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3. The evaluation of vA 
Now Ce,(€,q) is proportional to 
> (= 1A ag(q)ox(w), 
0 
where J,, as usual denotes Bessel functions of the first kind, w is vs cosh €, 
ind the coefficients are functions of g (= s/4). This expression is in turn 
proportional to ‘ 


r= Y Cy, (8)Jy,(w), (3.1) 


0 


where Cy, 1) D,,, and the D,, are the functions of s tabulated in 
2), pp. 41-57. It is clear from the tables that the C,, are positive for all 
positive s 

\s a preliminary step in the evaluation of a zero w, of (3.1) for a given 


ilue of s, equations (2.6) and (2.7) can be used to find a good approxima- 


on. Thus, since 1 A * for all eccentricities, we get wy .(w3—s) > a; 
it 1S ) . 9 ‘ 9 <« 
Wo st) \ {($8)* + x?! wi, Say, (3.2) 

where at = 33-4452. The fact that w)—w, is quite small (not greater 


in 0-18 for any s) is very helpful in the exact calculations of wy, particu 
ly for values of s greater than about unity. As the chosen value of s is 
reased from zero, it is found, during the process of evaluating w, exactly, 
increases from zero to a maximum of 0-176 when s is 25 and 

en decreases. When, for example, s is 100 (the largest value for which 
re tabulated), w, is 10-01 


ne coemecients in (2.5) 


and wy, is 10-136. In 
t, if it is assumed at this stage that (w5—s) > 7?/4 as s > 00 (see the end 


TNIS Se¢ 


ion), it follows quite easily from (3.2) that, for large values of s, 


9 > 
uw uw) ~~ 7°/SNS. 


‘eturning now to the precise determination of zeros of (3.1), it is clear 
hatfora positive valueofs lessthan unity, w,—w, isso small that theevalua- 
tion of the exact zero presents no difficulty; the function o is tabulated in 
e neighbourhood of the approximate zero w, and then w, is found by 
verse interpolation. For a larger value of s, the value of w,—w, found 


during the previous stage of the calculations serves as a useful estimate of 


the zero. Thus, when s = 16, w, = 4-227 and w,—, is a little greater than 
0-15, its value at s 13. The approximate values of o when w 4-39 and 


+40) are respectively 0-0038 and 0-0041, and they suggest tabulation 


round the value 4-395, as follows: 





“00002 209 "00059 747 *OOI75 170 
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Thus, when s = 16, w 4-39497 449 and, therefore. 


less sat 


toward 


f \8/Wp 0O-91013 042, 3-3 1580 O8. 


For certain suitable values of s ranging from zero to 100, the w-zeros of 


9» 


(3.1), and hence the pairs of corresponding values of e and Ay (1—e?), ar eRe 
evaluated in the manner described above. The results are contained j 5646 
Table 1. In each case o changes from positive to negative as w increases 49037 74 
through the zero; this is obvious when s is small since the function is they ‘aah 
dominated by the first term C,4,(w) and (, is slightly greater than unity sk 


The 



































TABLE J] mt 
piace ¢ 
| 9 ross 0 
0 \ o | yl e) | AVC e") 
" anianien ae SN — a = . 408 
674 
0°00 2°40452 | | | I 5°7931 240% 
0°02 2°40090 731 "99527 235 77320 Nur 
0°04 2°40899 425 | ‘99654 770 *763 : 
as ind in 
o’o 2°41315 360 "99310 747 74345 5 
4 
orl 11739 354 -98967 944 72379 15 lhese 
orl 42162 393 *95626 370 704260 
ee 42587 463 *98286 038 695450 77 
0°40 2°44742 S71 | *25541 632 | *96603 305 5°50990 73 
o*So0 |} 2°49195 300 | *35592 157 | *93336 701 $0997 
1°00 2°51494 409 | *39762 315 | °91754 882 194 3 
2°00 2°63568 177 | °5 "84355 924 1°94081 54 
,"00 2°70403 005 | *62063 950 | *77930 927 $°939 5 
4°00 2°89696 612 | *69037 742 | °72344 939 1°39241 27 
5°60 3°11 310 696 76015 118 64974 625 $°09143 49 
7°00 3°30110 576 | ‘80147 427 59802 926 )72 
10°00 3°05957 O57 «=| *S5701 551 | “51529 011 ,OISI47 
13°00 $°05430 973 | *88931 323 | °45729 857 3°43742 74 
16:00 1°39497 449 Q1013 042 | 41432 188 3°315800 
25°00 5°30150 516 13 . 7O 
30°00 24316 979 Of 2°97710 90 
50°00 27227 227 97 } 2°55594 4 
64°00 17493 080 9 | 49 3° 
100°00 I 640 1 | 74 
Y v : Y It ] 
x I 1/4 1740 II 
= — — he } 
_ : oe ‘ . . . , : \ | 
The value of A\/(1—e?) for s 100 is liable to an error of a unit or so! 


the fourth decimal place; this entry, however, is not necessary for subs 


quent calculations. Ip 
The values of A,(1—e?) at equal intervals in e, and in ,/(1—e?), are now WI 
constructed from Table I by the method of linear cross means to whic! 
reference has been made in section 2. (A method dependent upon divided ares 
differences is unsuitable owing to the very slow convergence of thes 
differences.) For example, in the calculation of the value of this functio! 1 
1 ine 


0-7 the convergence of the successive cross means proved to | 


when « 
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ENCIES OF 


ry than in all other cases; the convergence of the process 


finally assumed value of 4°35281 2 is shown in the Table LI. 


TABLE I] 


an 
88 | 4-3528 
ies are liabl be in error by two or three units in the sixth 
It is of interest to note here that the method of linear 
en ipplied yevaluate Ay (1 e7) when ¢ 1. gives the value 
sis correct to five figures for the theoretical limit 7*/4. 
lues of the principal frequency vA at one-tenth intervals in ¢ 
respectively are shown in columns 3 and 6 of Table II] 


ies can be given correct to four places of decimals only. 


TABLE II] 


\ ) Ay VA 
c | 4°2935 +7 
144 55 
"0% 17 
1507 C 7537 
I ( 7) 5°! 
int that over the major part of the range of the eccentricity 
frequent nereases remarkably slowly, and even when 
hat is whe! 0-995. vA is only 5-1309. 
evaluating when 0-1 yl e7) 0 is outlined below 
in ASW) ce formula 
{ Ow e (1, p. 389) 
Le,(E,,q) ~ ¢ | C_)COS Xo4 29 SIN Xo}; (3.o 
stant 
L/% sinh €,)/Z» tvs cosh7€,, 
Yo \ssinh €, — tan-"(tanh 3€,). 
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The zeros of (3.3) are such that 


cot x, = tan(d7— xp) —{1+-(1/a,)}-'!sinh €), 
that is 


Vs sinh €, = }7+-tan-!(tanh $£,)+ tan {sinh €,/(1-+-4vscosh7é,)!. (3.4 
to infinity, 2 +1 oe 
. Aj (1—e?) = ssinh?&, > 7?/4, 


as previously stated, and assumed in Table I. 


when 0-1 > ,/(1—e?) > 0. Thus, when ,/(1—e?) = 0-1, we get 


exp £) = 11/v99, 
and therefore 


vs sinh €, (8/99) ia+-tan-!(10— v99)+ tan-!{ v99/(99+ 400s)}, 


of this expression, are respectively: 
VS 15-62923 16°12755 16°14267 


VA 4°96729 5°12567 5:13048 


place. Thus the error incurred by use of the asymptotic formula over t! 


range of ,(1—e?) stated above is certain to be less than this amount. 


A note on a known formula for vA and its range of applicability 


in terms of the eccentricity: 


No greater claim is made for the formula than its accuracy for ‘smal 
eccentricities. However, using the expression (1—e?)! for the semi-min 
axis of an ellipse of area z, the following pairs of values can be derive 
from (3.5): 

é 0-4 0-6 O-8 

VA 2-40915 2°4331 2-5530 


These compare very favourably with the exact four-figure values in Tab 


in practical applications of the formula. That the formula can in fact! 
used for quite ‘large’ values of e is not surprising, now that the extreme 
slow rate of increase of VA over the major part of the full range of ¢! 
eccentricity has been established. 


Firstly, we note that the right side of (3.4) tends to 47. Hence, as s tend 


Secondly, we can apply (3.4) to evaluate vA with considerable accuracy 


approximately. By a method of successive approximations the consecutiv 


values of vs, and those of vA, obtained by using one, two, and three terms 


Comparison of the third approximation 5-1305 with the correspondin 


entry in Table II] shows an error of only four units in the fourth decima 


Herriot (6) has derived the following formula for this principal frequenc 


bvaA 2-4048... x (1 —0-25e2— 0-034640e4 — 0-010355e8 — 0-004650e8). (3.0 


III; even when e = 0-8, the error is probably small enough to be unimportat! 
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FREQ 


4, The evaluation of \, 
In this problem ¢d/en, and therefore ¢d/c&, are zero on the boundary 


nthe notation of section 2, v is equal to p, the principal eigenvalue, 

















| zeros are required of the function expressed in (2.9). 
; We use an expansion for ( e,(€,q) equivalent to that given in (1), p. 160; 
ve thee | le ised there are proportional, for a given s, to 
efhicie ( |)" D,,.., which are in effect tabulated in (2), 
aS—bDS ihe ¢ , are ositive for all positive Ss. Thus the function 
29) is DI ul ti 
Cul 
> (5) Js i vs cosh &,). 
ent the zero &, 0, which exists for all s, and applying 
ion for t! Bessel function derivative, it becomes cleat 
: we 1 juire zel f the expression 
ay 
SN ( I (uw) 2} Lar NS (wy) (4.1) 
range of values of s for which real zeros exist is, however, in this case 
ted | ondition ut \s/w, must not exceed unity, ie. that 
ondir 7 —§ 0. This condition restricts s to the range 0 < s < 3-56, 
eCcll eright s the inequality corresponding to ¢ i. So 0. Thus, after 
ver t ling (2.9) | und then making &, —- 0, we get the following equation 
nt - inn 
{) ‘ y°/12) 97q°/ 1152) 
ppre root is 0-890 and corresponds to the value 3-56 of s. 
quel \n alternat method pable of greater accuracy makes use of the 
* DD , s)cosh(2k je. 
: p rtional to C /). From this expression we derive our 
1-] . 
nf mely 
det 
iim Le “aa (2), 1)*De., 1 2. say. 
ppropriate value of s (i.e. 3-55928) is then found from the following 
f ot , \* 
n 1 
nt : : / 
I 4°5 5°0 
I %j cheat ees mee ee 
~ 14089 2 2 47618 4 
trel _ 
4 ° . 
of t ‘s Sincreases from zero over its limited range, the zero w, of (4.1) increases 
mp = 1-S411... to the slightly greater value ,/(3-55928) 1-8866. As 
Bb 
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tion and inverse interpolation, presents no special difficulty. Howevye 


within this range. Their values are given in the table below: 








| Cy Cs Cs 
3°51 1°13303 207 | 0°13815 87 | 0°00522 3 
3°52 | 348 799 865 71 257 
3°53 394 442 913 71 29 1 Cz = 000010 
54 | 440 134 962 73 32 5 
355 | 485 876 | o-rg011 80 359 








Column 2 of Table IV gives zeros of (4.1) for values of s from zero to 3:3; 
The other columns of the table give corresponding pairs of values of ; 


or of ,/(1—e?), and w3 = p/,/(1—e?), the intervals being, of course, unequa 


2 
2 


The values of w2 and e are liable to errors of 2 units in the seventh ar 
Th le f we and ‘e liable to e1 rz t tl th 
eighth places of decimals respectively.) 


TABLE IV 























I 2 ] 3 } 5 

= ‘'9 } te vu e*) to 
fohie) 184118 38 ° I 3°38995 78 
orl | 1°84229 59 | +17164 874 | -98515 822 | 3°39405 42 
0°25 1°84398 10 | ‘27115 247 | *96253 641 | 3°40026 59 
o°5 1°84683 40 | *38287 512 | ‘92380 O11 | 3°41079 58 
are) | 1°85271 21 | *53974 927 | *84182 583 | 3°43254 21 
2°0 1°86517 95 3°47559 46 
2°5 1°57175 12 52 37503560 49 
3°0 1°87863 50 | 92197 303 | °38725 409 | 3°52926 95 
3°5 | 1°88574 53 | °99208 O81 | *12553 O12 | 3°55603 53 
3°51 | 1°88589 of | -99342 979 | *11444 326 | 3°55658 15 
352 | 1°88603 83 | 99479 575 10218 172 | 3°55714 05 
3°53 | 1°88618 02 | -g9610 282 | -08819 957 | 3°55767 57 
3°54 1*88632 53 | °99743 601 | 07156 406 | 3°55822 31 
3°55 1'88647 05 | 99876 694 | 04964 474 | 355877 09 





The values of w2 and subsequent values of vu, at one-tenth interva! 
0 | I 


in e, and in ,/(1—e?), are here again evaluated by linear cross means, #! 
are shown in Table V. 

The nature of the convergence of the cross means, and differencing 
the final values, strongly suggest that the entries in columns 2 and 5 
Table V will have only slight errors in the fifth place of decimals, #! 
hence that the values of vj can confidently be given correct to four plac 
of decimals, as in the table. 


the latter interval is so small, the evaluation of wy, for a given s, by tabul. 


because the values of vu decrease quite rapidly to zero as s increases fron 
3:50 to its upper limit, sub-tabulation of the coefficients Cy, , is necessan 
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TABLE V 








3 | } 5 ¢ 
|| 2 2 | 

Vp jv e") 0 Ve 

| | 
5412 ] IO | 3*28996 I*S412 
309 o'9 3°41719 1°7537 
40 ° 3°44329 10507 
I¢ o'7 3°40753 | 1°5550 
o°O 3°49030 | 1°447! 


3°51029 | 1°3245 


o’4 1°1875 
7 0°3 1'0307 
o 0842 
orl 0°5964 
one) ° 





5. A note on the function (d/d£é)Ce,(€,¢) 
The function Ce,(€,q) involved in Problem I has a positive zero in € for 
ny assigned positive value of g, and consequently the principal frequency 
\inereases without limit as the eccentricity of the ellipse tends to unity. 
[he derivative of Ce,, however, has a positive zero in € only if 0 < q 0-89, 
nd thus \u is such that 8 Vu () for all values of the eccentricity. 
To illustrate the latter functional property, parametric curves of the 


tion (d/dé)Ce, are shown in Fig. 1, for s ty 0. 1. 3. 3. and 4. 


© 
n 


v 
mn 








rhe absolute values of this function depend on the way in which Ce, itself 


‘normalized; for convenience the normalization to be found in (2), p. XIX, 
s adopted hers Maes 


Ce, (0, q) .. 
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It then follows (pp. xx and xxii) that 


Ce,(€, 9) V(47)g, (s) 3 Cy 1 Jor. (\8 cosh &), 
0 
where y,, is the ‘joining factor’. Finally we have 
(d/dé)Ce, oa’, (4z)sinh € vsq,;, 
where o’ is the infinite series (4.1). 


The function vsg,, is tabulated in (2) (see p. 219). 


When s is zero, a’ is equal to } and vsg,, isin fact equal to 2, (2/7): her 
(d dé) e,(€, 0) sinh €, 
and this is in agreement with (2.9). 
The positive zeros of the function (5.1) for s 1, 2, and 3 are respectivel 


1-228, 0-779, and 0-409, Also it follows from (2.1) that the maxima oe 
for values of € satisfy ing 


(a bs cosh 2€)Ce, 0) ) 


i 2 _ , a 


where in this case a, is the ‘characteristic number’; the corresponding 
number as defined and tabulated in (2) (see p. 4) is be, a,+4s. Itis 

interest in connexion with (5.2) to note that the three appropriate valu 
of € are not zeros of Ce, but are found from the equation 


2a,/8 2(be,/s)—1. 


cosh 2E 
Fors 1.2,and3 we get, therefore, é 0-780, 0-466, and 0-238 respective 
at the maxima. The maximum values of the function could, if required, | 


found by direct substitution of s, €, and w vs cosh € /(be,) into (5.1 
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SUMMARY 


\ solut t 1 for t problem of an infinite beam under the action of 
vhere H he Heaviside unit function. The solution, which 

sed nko’s theor f transverse vibrations of beams, is obtained by 
se of | r transforms, fi which asymptotic approximations are obtained 
wtionary } Che solution shows that three wave-fronts aré 


uirection f1 tl loading point, and the stress distribution 


not differ greatly from that given by ordinary 


1. Notation 


dit t be section 


H(t) Heaviside unit function 
I m« ent of inertia of beam section 
is of gyration of beam section 


lass OT beam per unit length 


VW bending moment 
nve ion parametel ised in Fourier transforms 
Pr t insverse point oad 
q transverse distributed load 
time 
yy 
aistance Liong bean 
transverse deflexion of beam 
/ deflexion component due to bending 
/ deflexion component due to shear 
tr teeth 
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p density of beam material 
EI 
- ae 
T ct/k 
A wavelength 


2. Introduction 

The solution of the problem of an infinite uniform beam under the actior 
of a transverse point-load PH(t) has an important application in problems 
of transverse impact on finite beams. Solutions for similar infinite-bean 
problems have been obtained by Boussinesq (1) and others, using the 
ordinary equation for vibrations of beams, viz. 


Oty oy 


ox ot 


El 


This equation, however, leads to an incorrect solution, in which there is 


instantaneous propagation of disturbances along the beam, whereas it i: | 


known (Prescott (2)) that an exact solution to such a problem exhibits 
wave-fronts which are propagated with finite velocity. 

The solution given here by the author is based on Timoshenko’s mor 
accurate theory (3), (4), in which (1) is modified by the introduction o 
additional terms representing the effects of shear deflexion and rotar 
inertia. These modifications iead to a solution exhibiting wave-fronts 
similar to those expected in an exact solution. A similar problem, als 
using Timoshenko’s equation, has been investigated by Dengler a 
Goland (5) by means of Fourier—Laplace transforms. Their solutio1 
however, is in integral form, and not suitable for computation. In a 
earlier paper Robinson (6) gives a method for a step-by-step numerica 
solution which is, however, tedious except for investigating the initi 
spread of the disturbance. 

In the present paper the solution is obtained by means of Fourier trams 
forms, and, by using the method of stationary phase, asymptotic approx 
mations are obtained for the shear force and bending moment in the bea! 
These approximations have the advantage of giving directly the amplitud 
and wavelength of the stresses which, owing to dispersion, are hig! 
oscillatory. Special approximations are also obtained for the stresses 


the region of the wave-fronts. 


3. Equations of motion 
The equations of motion used are those derived by Timoshenko taking 


into account shear deflexion and rotary inertia, the total transver 
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leflexion being resolved into two components y, and y, due to bending 
snd shear respectively. 


The stress-strain relationships are 





M EIc% (2) 
Ua” 
for the bending moment and 
F os (3) 
Cu 


for the shear force, where C is the effective shear rigidity of the beam 
section, that is, the shear force required to produce unit shear strain at 
the neutral axis 

The equations of motion? of the beam are 


> 9 


‘on J. m =a (Yo +-Y.)- q(x) 0. (4) 
Ci" Ci" 
( 3) 3 
EI 20 _ mz 4 4 Cs — 9, (5) 
Cx? cauct= Ox 


\ simplification can be obtained by introducing non-dimensional 


riables, s, wu, and 7, defined by 


x 
Ss 
I: 
7 
u y | _ (6) 
k 
ct 
k 
Making these substitutions equations (2) to (5) become 
El C7u 
M rs (7) 
k cs 
Cu 
I ( (3) 
cs 
e2y 2 K3 
-(U, +Uu,) + — q(x) 0, (9) 
Os rT El 
ofu l cu 
b > of 8 (10) 
gv Cser* ‘ y OSs 
| - . . ° ° . 
the quantity EI/Ck* measures the effective ratio of flexural to shear 
tuliness. For a beam of rectangular section, y 4 approximately, but for 
In his f ierivation Timoshenko obtains a single equation for the total deflexion 


s) it mvenient, however, to treat the problem as one of coupled vibrations 
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flanged sections, which are designed for flexural rather than shear strength 
y is considerably larger. 
4. Fourier transforms 


Multiplying equations (9) and (10) by e’”§ and integrating with respect 


to s from —~» to +0, we get 
r 
n* - d2 , , 3 » Se 
u.-| = (WM), + Ug) - q(aye ps dg ) 1] 
974 dz? El 
— Oo 
a da or 
and pri, 4 2 —G, = 0, (12 
dr? vy 


/ 


' 
where @ is the Fourier transform of u defined by 


i(p, 7) | u(s, r)e?* ds. (13 
Also, because of the point load PH (t) applied at x — 0. we have it 
equation (11) ; 
re 2 
| q(x)e'’s ds - ] H(z). 14 
The complete solution of equations (11) and (12) satisfying the initia 
conditions @ 0 and di/dr Oatr 0 is 
: Pk? 
fy — pall + Aulpeosen t+ Aal pcos ever 
y Pke2 
ai. “___{] 1 B (p)cosw, 7+ B,(p)cos w,7}, 
E [p? ra ” 
where = 
A,(p) — 
Ws Wy 
w 
A,(p) —~- Se 
2 W5— Wy 
* yw 
B,(p) AE ee 3 
y(w5—w}j) 
yu 2 
Bp) = —1—£. 
2 y(w5 wy) J 
and w,, w, (ws > w,) are the positive roots of the equation 
(yw?— p*)(w* p")- w2 = 0. (1s 


The solution of the problem is now obtained by inversion of the Four’ 


transforms according to the usual inversion formula 


l ; 
u(s, 7) : a(p, rye"? dp. ;, 
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rengt! hus the bending moment and shear force are given by 
. Pk : dp 
VW | 1. (p)cos w, 7+ A,(p)cos w, tle~'8 f (20) 
on 2 2 p> 
P f dp 
| 5, (p)cos w,7 B,( p)cos wa, rh — 4 (21) 
Qat . . p 
The ce mn of the b is. of course, symmetrical about s 0 and 
her analysis will therefore be restricted to the region s 0. 
Iti e shown that, when s 7, the integrals (20) and (21) are identi- 


lhe integrands are single-valued and have no poles throughout 


plat id the pat f integration may therefore be deformed into 
circle in the lower half-plane. Along this path the integrals 
ll own to be ze vyhen s * 


5. Stationary-phase approximations—the general case 


Lient to resolve the bending moment into two components, viz. 
P ‘ dp 
1/ |-A,(p)cos w, TH : of (22 
| dp 
\/ | | P)COS Wy T é , 9 (23) 
p 
rie re 


I’ F dp 

F, | B,( p)cos w, 7} ps “I ’ (24) 
Zz Pp 
2 : d ) 

i | $,(p)cosws7e-'” la (25) 
2 J q " p 


‘\one of the integrands has a pole within the range of integration, and 
symptotic approximations for the integrals may therefore be obtained 


the method of stationary phase (Lamb (7)). 


msider the integral 


can dp 
i. pyexplif(p)| e 
‘ Pp 
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The approximation given by the stationary-phase method is isa VE 
A,(p,) - other 
5 1 l/s\ 7 l—$ aw . ‘| P 
V(2a) IF" (pa) 3 expli f(r) tz]. (28 
Pi 
The + or sign is taken according as f”(p,) is + or —, and p, is the 
value of p satisfying the stationary-phase condition 

f'(r,) = 0 iy | Equa 
range 


which may be written in the form 


dw s 
=-. (30) 
dp}, T 
The approximation (28) is valid sword’ that A,(p)/p? is not varying 
rapidly at p = p, and also that 
PP) SUP (Pv) (31 


{\@oy|\ 
Wap? Win 


Except for certain values “ p, this criterion is satisfied for large value 


; : | d3w 
that is, provided I 1 


dp* 


d*w, 











dp 


of + = ct/k and the solutions obtained are therefore asymptotic approxi 
mations, which are valid provided the disturbance has spread over : 
distance large compared with the depth of the beam. 

Using the procedure outlined above, the following approximations ar 





. TI 
obtained: 
uw Pk A (2 ») { d2 20, \\-3 - . . take’ 
m (277) “ | dp |} * COS(w, 7— ps 7) , (38 Figs 
NST | 
ies and 
Pk [A,(p){|d?w,| . j 
and M, = 2nr) aN 7 a * cos(w7— ps+ 17) . (34 cont 
oo } I P=Pp2 
where p, is the root of the equation It w 
a s 19° H 
goer — the 
dp Pp=Ppe ‘ P 

Similar approximations may be obtained for F, and F,, larit 
; a , ; and 
6. It will be understood that if the equation (30), for example, has tv - 
app 

real roots, then there will be two corresponding contributions, each of the 
form (33), to the approximation for /,. Should there be no real root, the 7. J 
the asymptotic approximation is zero. (' 


Fig. 1 shows the variation of dw,/dp and dw,/dp with p. These curves 
have been plotted for the particular case when y = 4, but the same essentia 
features are obtained in the general case when y > 1. It is seen that ther 
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TRANSIENT 


isa value of 8/7 (= a, say) for which equation (30) has repeated roots. For 


ther values of s/r the number of real roots of the equation are as follows: 


$/T x no root, 
\ sit >y two roots, 
sir >O > one root. 


Equation (35), on the other hand, has only one real root throughout the 








range ] $/7 0 
Lu, — 4 
ap 
~ du, 
‘i dp 
¥ f = — - a 
' i ie 
fy 
4 
l a . r 
2:0 p 
Fic. ] Variation of dw, dp, dw, dp, with Pp. 


The numerical evaluation of the above approximations has been under- 


ken for the particular case when y = 4, and the results are presented in 
Figs. 2 to 4, which show the amplitude and wavelength of M,,/(2ar)/Pk 
Wl,, (277)/ Pk plotted against s/r. The component M; is the additional 
ntribution due to the second root of equation (30) in the region 
x 8/T - Y 
it will be noted that the component M, is comparatively small. 

Here again, the same essential features of the solution are obtained in 
general case when 4 | and it is seen that there is a wave-front at 
ry * in addition to that at s — 7, together with an additional singu 

ity at s it. At all three points the condition (32) is not satisfied, 
| the general approximations hitherto obtained break down. Special 

pproximat will now be obtained for these points. 


/. Approximations when s ~ 7, s ~ try 


) 
Consider the integral 
dp 


5 (36) 
p° 


1,(p)expliw. r—ips 
a\) L[°Wo / 
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M, “Pk 
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207 
















0 T 
0-2 0-3 0 


» 
oO 
wn 


Fic. 2. Variation of amplitude of 7, /(2a7)/Pk with s/r ( v/ct). 
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curve gives the solution by ordinary beam theory.) 
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(The dotted 
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Fic. 3. Variation of amplitude of M, y(2a7)/Pk with s/t. 
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\\\ M. 
| 
~ 
» 
a 4 
\ a 
\ Rie 
™ le ay = 
ie. or 
es: 
M,’) y oo 
/ ue 
eth with s/r. (The dotted curve gives the 
dinary beam theory. 
, ©, and the bulk of the contribution to the integral occurs 
ro es of p. N when |p) is large, 
| 
, ‘ ) j 
2\p\(y—1) a 
(37) 
' l 
/ “~~ 
; l 
herefore, when the integral (36) has an asymptotic approxi- 
| iT dp 
Pp] epr Lps . -_ 4 (338) 
| 2p(9 1) Pp 
the use of the zero lower limit here is necessary to exclude an unwanted 
mtributi to the integra 38) in the region Pp 0. 
Similarly, when s > the bulk of the contribution to the integral 
‘ * dp 6 
1,(p)exp}| —tw.7 ips | . (39) 
2 p 
| = irs In the region p © and therefore an asymptotic approximation is 
07 dp 
ipr r (40) 
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Hence om 
2]. re . “ 

M, =~ pietie, | exp| ip s)+ =| ~ 

|! 2p(y—1)} p* 


4n(y—1) . 
=e LESA 
, s | 27(7—8) —— : 
F, = Pi lL MF =H 8). (43 


and similarly 





"F 
in the forms 


mm aU 


Orlevy—t)\" 
= Z| /(=5 —t 5) [Mle (45 


Equation (45) shows that there is a discontinuity of magnitude }P in th 





H(svy—r), (44 








shear force at s = ry-! 
8. Approximations when s \T 
When s = ar, equation (30) has repeated roots, p, = p; = &, say, al 
d*w,/dp* 0. 
Consider again the integral (26). Since both f’(é) and /"(€) are zero, w 
can write , 
f(p) = f(€)+= f’'(é)+... (46 


and therefore, when s = a7, an asymptotic approximation for the integra 
(26) is 


x 














é ° : eyarc¢ 
aN J expli fi é)| | exp UP 2 () dp 47 
Cc . fed 
- exp] Uf(ENSL'(E) | exp = | dfs, (48 





x 





thus it may be shown that, when s XT, 
Pk A 7|C3w, | 1 ( 
M, ~ 2-230 — 1?) ke - — | * cos(w, 7 ps) ; 4) 
2r| p? \2 | op | sal 
and similarly 
ys 
3 | “$INn(w, 7— ps) | ‘ y 


p=¢ 


Approximations for .W, and F; in the region s ~ ry~! may also be obtained 


The value of the integral in the expression (48) is 2-230 approximately, and 
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For the particular case when y = 4, a = 0-536, € = 0-844, and these 
approximations become 

M, ~ —1-26Pkr-* cos(0-3547—0-844s), (51) 


F, ~ —0-883Pr-* sin(0-3547— 0-844s). (52) 


9, Approximations when s is small 


Finally, when s/r > 0, p, — 0 and the approximation (33) breaks down 
wing to the presence of a pole at p = 0 in the integrand of (26). 
Now, when p is small, ° 
mm = Pp 
1 
A, a —| ;. (53) 
> ~ 
B, ] 


\symptotic approximations for V/, and F, in the region s,7 ~ Oare therefore 


s 


ret d 
M, ~ 5— | (1—cosptrje-ime— (54) 
27 p 
; Pf dies 
F, i ; (] COs p*r)e ps Py (55) 
mil p p 
The integra 55) is of standard form, giving 
: P : 7 
F, ~ —z{1—C(z)—S(z)}, (56) 
wh T s° a (57) 
where ; il 9 
4r 4ctk 
ind ¢ 2 S are Fresnel Ss integrals defined by 
rr 
. 4 (277) 


the integral (54) may be derived from (55) by integrating with respect to 


civin: 
Psk ; sin z — cosz 
u,= - I C(z)—S(z) + —_._—— \ (59) 
2 | J(2nz) | 
Expressions (56) and (59) are identical with those obtained for the shear 


lorce and bending moment by ordinary beam theory, using equation (1) 
Jones (8 
> ° ° . ’ ’ . . 
By using the isymptotic expansions for C(z), S(z) (Watson (9)), expression 


59) can be written in the following form, valid for large values of z: 


y  _2Ph (2) cos(e- |7). (60) 


\ 


(aT) 
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This asymptotic approximation is given by the dotted curves in Mig 
and 4. It is seen that, when y t, the solution by ordinary theory 
a good approximation (apart from a discrepancy in the wavelength) ing 


region 0 < s Ty 


/ 


10. Application to other loading conditions 

The formal extension of the present solution to problems in which# 
external load has any given time-variation may be achieved by theg 
of Duhamel’s integral. 

Of particular interest is the case when the beam is subjected to 
impulsive load. If the impulse is taken in the form of a Dirac 6-functigg 
the solution may be formally obtained by differentiating the preseq 
solution with respect to ¢. It follows that there will be an impulsive she 
force of magnitude }/, at the point s ry *, J, being the magnitude of 
applied impulse. From equations (41) and (44) it may also be shown tf 


there will be discontinuities in the bending moment of magnitude 


Iyevy 
ind ‘ 
2(y—1) 
From the point of view of technical applications, it is evident that, whe 
the beam is subjected to a sustained load, the maximum stresses W 


ultimately occur at the loading point and will be given satisfactorily § 


ordinary beam theory. Fora highly impulsive load, however, the maxim 


stresses are likely to occur at the wave fronts, in which case the ording 


theory will be inadequate. 
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